#include "custom_elements/stokes_3D_twofluid.h"

namespace Kratos {


void Stokes3DTwoFluid::ComputeGaussPointLHSContribution(BoundedMatrix<double,16,16>& lhs, const element_data<4,3>& data)
    {
        const int nnodes = 4;
        const int dim = 3;

        const double rho = inner_prod(data.N, data.rho);
        const double& bdf0 = data.bdf0;

        //get constitutive matrix
        const Matrix& C = data.C;

        //get shape function values
        const BoundedMatrix<double,nnodes,dim>& DN = data.DN_DX;
        const array_1d<double,nnodes>& N = data.N;


        //compute an equivalent tau by Bitrans*c*Bi
        const double tau_denom = (C(3,3) + C(4,4) + C(5,5))*(pow(DN(0,0), 2) + pow(DN(0,1), 2) + pow(DN(0,2), 2) + pow(DN(1,0), 2) + pow(DN(1,1), 2) + pow(DN(1,2), 2) + pow(DN(2,0), 2) + pow(DN(2,1), 2) + pow(DN(2,2), 2) + pow(DN(3,0), 2) + pow(DN(3,1), 2) + pow(DN(3,2), 2));

        const double tau1 = 1.0/(tau_denom*rho);
        const double tau2 = (C(3,3) + C(4,4) + C(5,5))/(6.0*rho);

        const double DN_0_0 = DN(0,0); const double DN_0_1 = DN(0,1); const double DN_0_2 = DN(0,2);
        const double DN_1_0 = DN(1,0); const double DN_1_1 = DN(1,1); const double DN_1_2 = DN(1,2);
        const double DN_2_0 = DN(2,0); const double DN_2_1 = DN(2,1); const double DN_2_2 = DN(2,2);
        const double DN_3_0 = DN(3,0); const double DN_3_1 = DN(3,1); const double DN_3_2 = DN(3,2);

        const double N_0 = N[0];
        const double N_1 = N[1];
        const double N_2 = N[2];
        const double N_3 = N[3];

        const double C_0_0 = C(0,0); const double C_0_1 = C(0,1); const double C_0_2 = C(0,2); const double C_0_3 = C(0,3); const double C_0_4 = C(0,4); const double C_0_5 = C(0,5);
        const double C_1_1 = C(1,1); const double C_1_2 = C(1,2); const double C_1_3 = C(1,3); const double C_1_4 = C(1,4); const double C_1_5 = C(1,5);
        const double C_2_2 = C(2,2); const double C_2_3 = C(2,3); const double C_2_4 = C(2,4); const double C_2_5 = C(2,5);
        const double C_3_3 = C(3,3); const double C_3_4 = C(3,4); const double C_3_5 = C(3,5);
        const double C_4_4 = C(4,4); const double C_4_5 = C(4,5); const double C_5_5 = C(5,5);

        const double clhs0 = bdf0*rho;
const double clhs1 = pow(N_0, 2)*clhs0;
const double clhs2 = pow(DN_0_0, 2);
const double clhs3 = rho*tau2;
const double clhs4 = C_0_0*DN_0_0 + C_0_3*DN_0_1 + C_0_5*DN_0_2;
const double clhs5 = DN_0_0*clhs4;
const double clhs6 = C_0_3*DN_0_0;
const double clhs7 = C_3_3*DN_0_1 + C_3_5*DN_0_2 + clhs6;
const double clhs8 = DN_0_1*clhs7;
const double clhs9 = C_0_5*DN_0_0;
const double clhs10 = C_3_5*DN_0_1 + C_5_5*DN_0_2 + clhs9;
const double clhs11 = DN_0_2*clhs10;
const double clhs12 = clhs0*tau1;
const double clhs13 = DN_0_0*clhs3;
const double clhs14 = -DN_0_1*clhs13;
const double clhs15 = C_0_1*DN_0_1 + C_0_4*DN_0_2 + clhs6;
const double clhs16 = DN_0_0*clhs15;
const double clhs17 = C_1_3*DN_0_1;
const double clhs18 = C_3_3*DN_0_0 + C_3_4*DN_0_2 + clhs17;
const double clhs19 = DN_0_1*clhs18;
const double clhs20 = C_3_5*DN_0_0;
const double clhs21 = C_4_5*DN_0_2;
const double clhs22 = C_1_5*DN_0_1 + clhs20 + clhs21;
const double clhs23 = DN_0_2*clhs22;
const double clhs24 = DN_0_0*clhs7;
const double clhs25 = C_0_1*DN_0_0 + C_1_5*DN_0_2 + clhs17;
const double clhs26 = DN_0_1*clhs25;
const double clhs27 = C_3_4*DN_0_1;
const double clhs28 = C_0_4*DN_0_0 + clhs21 + clhs27;
const double clhs29 = DN_0_2*clhs28;
const double clhs30 = -DN_0_2*clhs13;
const double clhs31 = C_0_2*DN_0_2 + C_0_4*DN_0_1 + clhs9;
const double clhs32 = DN_0_0*clhs31;
const double clhs33 = C_2_3*DN_0_2 + clhs20 + clhs27;
const double clhs34 = DN_0_1*clhs33;
const double clhs35 = C_2_5*DN_0_2;
const double clhs36 = C_4_5*DN_0_1 + C_5_5*DN_0_0 + clhs35;
const double clhs37 = DN_0_2*clhs36;
const double clhs38 = DN_0_0*clhs10;
const double clhs39 = DN_0_1*clhs28;
const double clhs40 = C_0_2*DN_0_0 + C_2_3*DN_0_1 + clhs35;
const double clhs41 = DN_0_2*clhs40;
const double clhs42 = DN_0_0*N_0;
const double clhs43 = C_0_0*DN_1_0 + C_0_3*DN_1_1 + C_0_5*DN_1_2;
const double clhs44 = DN_0_0*clhs43;
const double clhs45 = C_0_3*DN_1_0;
const double clhs46 = C_3_3*DN_1_1 + C_3_5*DN_1_2 + clhs45;
const double clhs47 = DN_0_1*clhs46;
const double clhs48 = C_0_5*DN_1_0;
const double clhs49 = C_3_5*DN_1_1 + C_5_5*DN_1_2 + clhs48;
const double clhs50 = DN_0_2*clhs49;
const double clhs51 = DN_1_0*clhs4;
const double clhs52 = DN_1_1*clhs7;
const double clhs53 = DN_1_2*clhs10;
const double clhs54 = N_0*clhs0;
const double clhs55 = N_1*clhs54;
const double clhs56 = DN_0_0*DN_1_0;
const double clhs57 = -clhs3*clhs56 + clhs55;
const double clhs58 = -DN_1_1*clhs13;
const double clhs59 = C_0_1*DN_1_1 + C_0_4*DN_1_2 + clhs45;
const double clhs60 = DN_0_0*clhs59;
const double clhs61 = C_1_3*DN_1_1;
const double clhs62 = C_3_3*DN_1_0 + C_3_4*DN_1_2 + clhs61;
const double clhs63 = DN_0_1*clhs62;
const double clhs64 = C_3_5*DN_1_0;
const double clhs65 = C_4_5*DN_1_2;
const double clhs66 = C_1_5*DN_1_1 + clhs64 + clhs65;
const double clhs67 = DN_0_2*clhs66;
const double clhs68 = DN_1_0*clhs7;
const double clhs69 = DN_1_1*clhs25;
const double clhs70 = DN_1_2*clhs28;
const double clhs71 = -DN_1_2*clhs13;
const double clhs72 = C_0_2*DN_1_2 + C_0_4*DN_1_1 + clhs48;
const double clhs73 = DN_0_0*clhs72;
const double clhs74 = C_3_4*DN_1_1;
const double clhs75 = C_2_3*DN_1_2 + clhs64 + clhs74;
const double clhs76 = DN_0_1*clhs75;
const double clhs77 = C_2_5*DN_1_2;
const double clhs78 = C_4_5*DN_1_1 + C_5_5*DN_1_0 + clhs77;
const double clhs79 = DN_0_2*clhs78;
const double clhs80 = DN_1_0*clhs10;
const double clhs81 = DN_1_1*clhs28;
const double clhs82 = DN_1_2*clhs40;
const double clhs83 = DN_0_0*N_1;
const double clhs84 = C_0_0*DN_2_0 + C_0_3*DN_2_1 + C_0_5*DN_2_2;
const double clhs85 = DN_0_0*clhs84;
const double clhs86 = C_0_3*DN_2_0;
const double clhs87 = C_3_3*DN_2_1 + C_3_5*DN_2_2 + clhs86;
const double clhs88 = DN_0_1*clhs87;
const double clhs89 = C_0_5*DN_2_0;
const double clhs90 = C_3_5*DN_2_1 + C_5_5*DN_2_2 + clhs89;
const double clhs91 = DN_0_2*clhs90;
const double clhs92 = DN_2_0*clhs4;
const double clhs93 = DN_2_1*clhs7;
const double clhs94 = DN_2_2*clhs10;
const double clhs95 = N_2*clhs54;
const double clhs96 = DN_0_0*DN_2_0;
const double clhs97 = -clhs3*clhs96 + clhs95;
const double clhs98 = -DN_2_1*clhs13;
const double clhs99 = C_0_1*DN_2_1 + C_0_4*DN_2_2 + clhs86;
const double clhs100 = DN_0_0*clhs99;
const double clhs101 = C_1_3*DN_2_1;
const double clhs102 = C_3_3*DN_2_0 + C_3_4*DN_2_2 + clhs101;
const double clhs103 = DN_0_1*clhs102;
const double clhs104 = C_3_5*DN_2_0;
const double clhs105 = C_4_5*DN_2_2;
const double clhs106 = C_1_5*DN_2_1 + clhs104 + clhs105;
const double clhs107 = DN_0_2*clhs106;
const double clhs108 = DN_2_0*clhs7;
const double clhs109 = DN_2_1*clhs25;
const double clhs110 = DN_2_2*clhs28;
const double clhs111 = -DN_2_2*clhs13;
const double clhs112 = C_0_2*DN_2_2 + C_0_4*DN_2_1 + clhs89;
const double clhs113 = DN_0_0*clhs112;
const double clhs114 = C_3_4*DN_2_1;
const double clhs115 = C_2_3*DN_2_2 + clhs104 + clhs114;
const double clhs116 = DN_0_1*clhs115;
const double clhs117 = C_2_5*DN_2_2;
const double clhs118 = C_4_5*DN_2_1 + C_5_5*DN_2_0 + clhs117;
const double clhs119 = DN_0_2*clhs118;
const double clhs120 = DN_2_0*clhs10;
const double clhs121 = DN_2_1*clhs28;
const double clhs122 = DN_2_2*clhs40;
const double clhs123 = DN_0_0*N_2;
const double clhs124 = C_0_0*DN_3_0 + C_0_3*DN_3_1 + C_0_5*DN_3_2;
const double clhs125 = DN_0_0*clhs124;
const double clhs126 = C_0_3*DN_3_0;
const double clhs127 = C_3_3*DN_3_1 + C_3_5*DN_3_2 + clhs126;
const double clhs128 = DN_0_1*clhs127;
const double clhs129 = C_0_5*DN_3_0;
const double clhs130 = C_3_5*DN_3_1 + C_5_5*DN_3_2 + clhs129;
const double clhs131 = DN_0_2*clhs130;
const double clhs132 = DN_3_0*clhs4;
const double clhs133 = DN_3_1*clhs7;
const double clhs134 = DN_3_2*clhs10;
const double clhs135 = N_3*clhs54;
const double clhs136 = DN_0_0*DN_3_0;
const double clhs137 = clhs135 - clhs136*clhs3;
const double clhs138 = -DN_3_1*clhs13;
const double clhs139 = C_0_1*DN_3_1 + C_0_4*DN_3_2 + clhs126;
const double clhs140 = DN_0_0*clhs139;
const double clhs141 = C_1_3*DN_3_1;
const double clhs142 = C_3_3*DN_3_0 + C_3_4*DN_3_2 + clhs141;
const double clhs143 = DN_0_1*clhs142;
const double clhs144 = C_3_5*DN_3_0;
const double clhs145 = C_4_5*DN_3_2;
const double clhs146 = C_1_5*DN_3_1 + clhs144 + clhs145;
const double clhs147 = DN_0_2*clhs146;
const double clhs148 = DN_3_0*clhs7;
const double clhs149 = DN_3_1*clhs25;
const double clhs150 = DN_3_2*clhs28;
const double clhs151 = -DN_3_2*clhs13;
const double clhs152 = C_0_2*DN_3_2 + C_0_4*DN_3_1 + clhs129;
const double clhs153 = DN_0_0*clhs152;
const double clhs154 = C_3_4*DN_3_1;
const double clhs155 = C_2_3*DN_3_2 + clhs144 + clhs154;
const double clhs156 = DN_0_1*clhs155;
const double clhs157 = C_2_5*DN_3_2;
const double clhs158 = C_4_5*DN_3_1 + C_5_5*DN_3_0 + clhs157;
const double clhs159 = DN_0_2*clhs158;
const double clhs160 = DN_3_0*clhs10;
const double clhs161 = DN_3_1*clhs28;
const double clhs162 = DN_3_2*clhs40;
const double clhs163 = DN_0_0*N_3;
const double clhs164 = pow(DN_0_1, 2);
const double clhs165 = DN_0_0*clhs18;
const double clhs166 = C_1_1*DN_0_1 + C_1_3*DN_0_0 + C_1_4*DN_0_2;
const double clhs167 = DN_0_1*clhs166;
const double clhs168 = C_1_4*DN_0_1;
const double clhs169 = C_3_4*DN_0_0 + C_4_4*DN_0_2 + clhs168;
const double clhs170 = DN_0_2*clhs169;
const double clhs171 = DN_0_1*clhs3;
const double clhs172 = -DN_0_2*clhs171;
const double clhs173 = DN_0_0*clhs33;
const double clhs174 = C_1_2*DN_0_2 + C_1_5*DN_0_0 + clhs168;
const double clhs175 = DN_0_1*clhs174;
const double clhs176 = C_2_4*DN_0_2;
const double clhs177 = C_4_4*DN_0_1 + C_4_5*DN_0_0 + clhs176;
const double clhs178 = DN_0_2*clhs177;
const double clhs179 = DN_0_0*clhs22;
const double clhs180 = DN_0_1*clhs169;
const double clhs181 = C_1_2*DN_0_1 + C_2_3*DN_0_0 + clhs176;
const double clhs182 = DN_0_2*clhs181;
const double clhs183 = DN_0_1*N_0;
const double clhs184 = -DN_1_0*clhs171;
const double clhs185 = DN_0_0*clhs46;
const double clhs186 = C_0_1*DN_1_0 + C_1_5*DN_1_2 + clhs61;
const double clhs187 = DN_0_1*clhs186;
const double clhs188 = C_0_4*DN_1_0 + clhs65 + clhs74;
const double clhs189 = DN_0_2*clhs188;
const double clhs190 = DN_1_0*clhs15;
const double clhs191 = DN_1_1*clhs18;
const double clhs192 = DN_1_2*clhs22;
const double clhs193 = DN_0_0*clhs62;
const double clhs194 = C_1_1*DN_1_1 + C_1_3*DN_1_0 + C_1_4*DN_1_2;
const double clhs195 = DN_0_1*clhs194;
const double clhs196 = C_1_4*DN_1_1;
const double clhs197 = C_3_4*DN_1_0 + C_4_4*DN_1_2 + clhs196;
const double clhs198 = DN_0_2*clhs197;
const double clhs199 = DN_1_0*clhs18;
const double clhs200 = DN_1_1*clhs166;
const double clhs201 = DN_1_2*clhs169;
const double clhs202 = DN_0_1*DN_1_1;
const double clhs203 = -clhs202*clhs3 + clhs55;
const double clhs204 = -DN_1_2*clhs171;
const double clhs205 = DN_0_0*clhs75;
const double clhs206 = C_1_2*DN_1_2 + C_1_5*DN_1_0 + clhs196;
const double clhs207 = DN_0_1*clhs206;
const double clhs208 = C_2_4*DN_1_2;
const double clhs209 = C_4_4*DN_1_1 + C_4_5*DN_1_0 + clhs208;
const double clhs210 = DN_0_2*clhs209;
const double clhs211 = DN_1_0*clhs22;
const double clhs212 = DN_1_1*clhs169;
const double clhs213 = DN_1_2*clhs181;
const double clhs214 = DN_0_1*N_1;
const double clhs215 = -DN_2_0*clhs171;
const double clhs216 = DN_0_0*clhs87;
const double clhs217 = C_0_1*DN_2_0 + C_1_5*DN_2_2 + clhs101;
const double clhs218 = DN_0_1*clhs217;
const double clhs219 = C_0_4*DN_2_0 + clhs105 + clhs114;
const double clhs220 = DN_0_2*clhs219;
const double clhs221 = DN_2_0*clhs15;
const double clhs222 = DN_2_1*clhs18;
const double clhs223 = DN_2_2*clhs22;
const double clhs224 = DN_0_0*clhs102;
const double clhs225 = C_1_1*DN_2_1 + C_1_3*DN_2_0 + C_1_4*DN_2_2;
const double clhs226 = DN_0_1*clhs225;
const double clhs227 = C_1_4*DN_2_1;
const double clhs228 = C_3_4*DN_2_0 + C_4_4*DN_2_2 + clhs227;
const double clhs229 = DN_0_2*clhs228;
const double clhs230 = DN_2_0*clhs18;
const double clhs231 = DN_2_1*clhs166;
const double clhs232 = DN_2_2*clhs169;
const double clhs233 = DN_0_1*DN_2_1;
const double clhs234 = -clhs233*clhs3 + clhs95;
const double clhs235 = -DN_2_2*clhs171;
const double clhs236 = DN_0_0*clhs115;
const double clhs237 = C_1_2*DN_2_2 + C_1_5*DN_2_0 + clhs227;
const double clhs238 = DN_0_1*clhs237;
const double clhs239 = C_2_4*DN_2_2;
const double clhs240 = C_4_4*DN_2_1 + C_4_5*DN_2_0 + clhs239;
const double clhs241 = DN_0_2*clhs240;
const double clhs242 = DN_2_0*clhs22;
const double clhs243 = DN_2_1*clhs169;
const double clhs244 = DN_2_2*clhs181;
const double clhs245 = DN_0_1*N_2;
const double clhs246 = -DN_3_0*clhs171;
const double clhs247 = DN_0_0*clhs127;
const double clhs248 = C_0_1*DN_3_0 + C_1_5*DN_3_2 + clhs141;
const double clhs249 = DN_0_1*clhs248;
const double clhs250 = C_0_4*DN_3_0 + clhs145 + clhs154;
const double clhs251 = DN_0_2*clhs250;
const double clhs252 = DN_3_0*clhs15;
const double clhs253 = DN_3_1*clhs18;
const double clhs254 = DN_3_2*clhs22;
const double clhs255 = DN_0_0*clhs142;
const double clhs256 = C_1_1*DN_3_1 + C_1_3*DN_3_0 + C_1_4*DN_3_2;
const double clhs257 = DN_0_1*clhs256;
const double clhs258 = C_1_4*DN_3_1;
const double clhs259 = C_3_4*DN_3_0 + C_4_4*DN_3_2 + clhs258;
const double clhs260 = DN_0_2*clhs259;
const double clhs261 = DN_3_0*clhs18;
const double clhs262 = DN_3_1*clhs166;
const double clhs263 = DN_3_2*clhs169;
const double clhs264 = DN_0_1*DN_3_1;
const double clhs265 = clhs135 - clhs264*clhs3;
const double clhs266 = -DN_3_2*clhs171;
const double clhs267 = DN_0_0*clhs155;
const double clhs268 = C_1_2*DN_3_2 + C_1_5*DN_3_0 + clhs258;
const double clhs269 = DN_0_1*clhs268;
const double clhs270 = C_2_4*DN_3_2;
const double clhs271 = C_4_4*DN_3_1 + C_4_5*DN_3_0 + clhs270;
const double clhs272 = DN_0_2*clhs271;
const double clhs273 = DN_3_0*clhs22;
const double clhs274 = DN_3_1*clhs169;
const double clhs275 = DN_3_2*clhs181;
const double clhs276 = DN_0_1*N_3;
const double clhs277 = pow(DN_0_2, 2);
const double clhs278 = DN_0_0*clhs36;
const double clhs279 = DN_0_1*clhs177;
const double clhs280 = C_2_2*DN_0_2 + C_2_4*DN_0_1 + C_2_5*DN_0_0;
const double clhs281 = DN_0_2*clhs280;
const double clhs282 = DN_0_2*N_0;
const double clhs283 = DN_0_2*clhs3;
const double clhs284 = -DN_1_0*clhs283;
const double clhs285 = DN_0_0*clhs49;
const double clhs286 = DN_0_1*clhs188;
const double clhs287 = C_0_2*DN_1_0 + C_2_3*DN_1_1 + clhs77;
const double clhs288 = DN_0_2*clhs287;
const double clhs289 = DN_1_0*clhs31;
const double clhs290 = DN_1_1*clhs33;
const double clhs291 = DN_1_2*clhs36;
const double clhs292 = -DN_1_1*clhs283;
const double clhs293 = DN_0_0*clhs66;
const double clhs294 = DN_0_1*clhs197;
const double clhs295 = C_1_2*DN_1_1 + C_2_3*DN_1_0 + clhs208;
const double clhs296 = DN_0_2*clhs295;
const double clhs297 = DN_1_0*clhs33;
const double clhs298 = DN_1_1*clhs174;
const double clhs299 = DN_1_2*clhs177;
const double clhs300 = DN_0_0*clhs78;
const double clhs301 = DN_0_1*clhs209;
const double clhs302 = C_2_2*DN_1_2 + C_2_4*DN_1_1 + C_2_5*DN_1_0;
const double clhs303 = DN_0_2*clhs302;
const double clhs304 = DN_1_0*clhs36;
const double clhs305 = DN_1_1*clhs177;
const double clhs306 = DN_1_2*clhs280;
const double clhs307 = DN_0_2*DN_1_2;
const double clhs308 = -clhs3*clhs307 + clhs55;
const double clhs309 = DN_0_2*N_1;
const double clhs310 = -DN_2_0*clhs283;
const double clhs311 = DN_0_0*clhs90;
const double clhs312 = DN_0_1*clhs219;
const double clhs313 = C_0_2*DN_2_0 + C_2_3*DN_2_1 + clhs117;
const double clhs314 = DN_0_2*clhs313;
const double clhs315 = DN_2_0*clhs31;
const double clhs316 = DN_2_1*clhs33;
const double clhs317 = DN_2_2*clhs36;
const double clhs318 = -DN_2_1*clhs283;
const double clhs319 = DN_0_0*clhs106;
const double clhs320 = DN_0_1*clhs228;
const double clhs321 = C_1_2*DN_2_1 + C_2_3*DN_2_0 + clhs239;
const double clhs322 = DN_0_2*clhs321;
const double clhs323 = DN_2_0*clhs33;
const double clhs324 = DN_2_1*clhs174;
const double clhs325 = DN_2_2*clhs177;
const double clhs326 = DN_0_0*clhs118;
const double clhs327 = DN_0_1*clhs240;
const double clhs328 = C_2_2*DN_2_2 + C_2_4*DN_2_1 + C_2_5*DN_2_0;
const double clhs329 = DN_0_2*clhs328;
const double clhs330 = DN_2_0*clhs36;
const double clhs331 = DN_2_1*clhs177;
const double clhs332 = DN_2_2*clhs280;
const double clhs333 = DN_0_2*DN_2_2;
const double clhs334 = -clhs3*clhs333 + clhs95;
const double clhs335 = DN_0_2*N_2;
const double clhs336 = -DN_3_0*clhs283;
const double clhs337 = DN_0_0*clhs130;
const double clhs338 = DN_0_1*clhs250;
const double clhs339 = C_0_2*DN_3_0 + C_2_3*DN_3_1 + clhs157;
const double clhs340 = DN_0_2*clhs339;
const double clhs341 = DN_3_0*clhs31;
const double clhs342 = DN_3_1*clhs33;
const double clhs343 = DN_3_2*clhs36;
const double clhs344 = -DN_3_1*clhs283;
const double clhs345 = DN_0_0*clhs146;
const double clhs346 = DN_0_1*clhs259;
const double clhs347 = C_1_2*DN_3_1 + C_2_3*DN_3_0 + clhs270;
const double clhs348 = DN_0_2*clhs347;
const double clhs349 = DN_3_0*clhs33;
const double clhs350 = DN_3_1*clhs174;
const double clhs351 = DN_3_2*clhs177;
const double clhs352 = DN_0_0*clhs158;
const double clhs353 = DN_0_1*clhs271;
const double clhs354 = C_2_2*DN_3_2 + C_2_4*DN_3_1 + C_2_5*DN_3_0;
const double clhs355 = DN_0_2*clhs354;
const double clhs356 = DN_3_0*clhs36;
const double clhs357 = DN_3_1*clhs177;
const double clhs358 = DN_3_2*clhs280;
const double clhs359 = DN_0_2*DN_3_2;
const double clhs360 = clhs135 - clhs3*clhs359;
const double clhs361 = DN_0_2*N_3;
const double clhs362 = bdf0*pow(rho, 2)*tau1;
const double clhs363 = clhs362 + 1;
const double clhs364 = rho*tau1;
const double clhs365 = DN_1_0*N_0;
const double clhs366 = DN_1_1*N_0;
const double clhs367 = DN_1_2*N_0;
const double clhs368 = clhs364*(clhs202 + clhs307 + clhs56);
const double clhs369 = DN_2_0*N_0;
const double clhs370 = DN_2_1*N_0;
const double clhs371 = DN_2_2*N_0;
const double clhs372 = clhs364*(clhs233 + clhs333 + clhs96);
const double clhs373 = DN_3_0*N_0;
const double clhs374 = DN_3_1*N_0;
const double clhs375 = DN_3_2*N_0;
const double clhs376 = clhs364*(clhs136 + clhs264 + clhs359);
const double clhs377 = pow(N_1, 2)*clhs0;
const double clhs378 = pow(DN_1_0, 2);
const double clhs379 = DN_1_0*clhs43;
const double clhs380 = DN_1_1*clhs46;
const double clhs381 = DN_1_2*clhs49;
const double clhs382 = DN_1_0*clhs3;
const double clhs383 = -DN_1_1*clhs382;
const double clhs384 = DN_1_0*clhs59;
const double clhs385 = DN_1_1*clhs62;
const double clhs386 = DN_1_2*clhs66;
const double clhs387 = DN_1_0*clhs46;
const double clhs388 = DN_1_1*clhs186;
const double clhs389 = DN_1_2*clhs188;
const double clhs390 = -DN_1_2*clhs382;
const double clhs391 = DN_1_0*clhs72;
const double clhs392 = DN_1_1*clhs75;
const double clhs393 = DN_1_2*clhs78;
const double clhs394 = DN_1_0*clhs49;
const double clhs395 = DN_1_1*clhs188;
const double clhs396 = DN_1_2*clhs287;
const double clhs397 = DN_1_0*N_1;
const double clhs398 = DN_1_0*clhs84;
const double clhs399 = DN_1_1*clhs87;
const double clhs400 = DN_1_2*clhs90;
const double clhs401 = DN_2_0*clhs43;
const double clhs402 = DN_2_1*clhs46;
const double clhs403 = DN_2_2*clhs49;
const double clhs404 = N_1*clhs0;
const double clhs405 = N_2*clhs404;
const double clhs406 = DN_1_0*DN_2_0;
const double clhs407 = -clhs3*clhs406 + clhs405;
const double clhs408 = -DN_2_1*clhs382;
const double clhs409 = DN_1_0*clhs99;
const double clhs410 = DN_1_1*clhs102;
const double clhs411 = DN_1_2*clhs106;
const double clhs412 = DN_2_0*clhs46;
const double clhs413 = DN_2_1*clhs186;
const double clhs414 = DN_2_2*clhs188;
const double clhs415 = -DN_2_2*clhs382;
const double clhs416 = DN_1_0*clhs112;
const double clhs417 = DN_1_1*clhs115;
const double clhs418 = DN_1_2*clhs118;
const double clhs419 = DN_2_0*clhs49;
const double clhs420 = DN_2_1*clhs188;
const double clhs421 = DN_2_2*clhs287;
const double clhs422 = DN_1_0*N_2;
const double clhs423 = DN_1_0*clhs124;
const double clhs424 = DN_1_1*clhs127;
const double clhs425 = DN_1_2*clhs130;
const double clhs426 = DN_3_0*clhs43;
const double clhs427 = DN_3_1*clhs46;
const double clhs428 = DN_3_2*clhs49;
const double clhs429 = N_3*clhs404;
const double clhs430 = DN_1_0*DN_3_0;
const double clhs431 = -clhs3*clhs430 + clhs429;
const double clhs432 = -DN_3_1*clhs382;
const double clhs433 = DN_1_0*clhs139;
const double clhs434 = DN_1_1*clhs142;
const double clhs435 = DN_1_2*clhs146;
const double clhs436 = DN_3_0*clhs46;
const double clhs437 = DN_3_1*clhs186;
const double clhs438 = DN_3_2*clhs188;
const double clhs439 = -DN_3_2*clhs382;
const double clhs440 = DN_1_0*clhs152;
const double clhs441 = DN_1_1*clhs155;
const double clhs442 = DN_1_2*clhs158;
const double clhs443 = DN_3_0*clhs49;
const double clhs444 = DN_3_1*clhs188;
const double clhs445 = DN_3_2*clhs287;
const double clhs446 = DN_1_0*N_3;
const double clhs447 = pow(DN_1_1, 2);
const double clhs448 = DN_1_0*clhs62;
const double clhs449 = DN_1_1*clhs194;
const double clhs450 = DN_1_2*clhs197;
const double clhs451 = DN_1_1*clhs3;
const double clhs452 = -DN_1_2*clhs451;
const double clhs453 = DN_1_0*clhs75;
const double clhs454 = DN_1_1*clhs206;
const double clhs455 = DN_1_2*clhs209;
const double clhs456 = DN_1_0*clhs66;
const double clhs457 = DN_1_1*clhs197;
const double clhs458 = DN_1_2*clhs295;
const double clhs459 = DN_1_1*N_1;
const double clhs460 = -DN_2_0*clhs451;
const double clhs461 = DN_1_0*clhs87;
const double clhs462 = DN_1_1*clhs217;
const double clhs463 = DN_1_2*clhs219;
const double clhs464 = DN_2_0*clhs59;
const double clhs465 = DN_2_1*clhs62;
const double clhs466 = DN_2_2*clhs66;
const double clhs467 = DN_1_0*clhs102;
const double clhs468 = DN_1_1*clhs225;
const double clhs469 = DN_1_2*clhs228;
const double clhs470 = DN_2_0*clhs62;
const double clhs471 = DN_2_1*clhs194;
const double clhs472 = DN_2_2*clhs197;
const double clhs473 = DN_1_1*DN_2_1;
const double clhs474 = -clhs3*clhs473 + clhs405;
const double clhs475 = -DN_2_2*clhs451;
const double clhs476 = DN_1_0*clhs115;
const double clhs477 = DN_1_1*clhs237;
const double clhs478 = DN_1_2*clhs240;
const double clhs479 = DN_2_0*clhs66;
const double clhs480 = DN_2_1*clhs197;
const double clhs481 = DN_2_2*clhs295;
const double clhs482 = DN_1_1*N_2;
const double clhs483 = -DN_3_0*clhs451;
const double clhs484 = DN_1_0*clhs127;
const double clhs485 = DN_1_1*clhs248;
const double clhs486 = DN_1_2*clhs250;
const double clhs487 = DN_3_0*clhs59;
const double clhs488 = DN_3_1*clhs62;
const double clhs489 = DN_3_2*clhs66;
const double clhs490 = DN_1_0*clhs142;
const double clhs491 = DN_1_1*clhs256;
const double clhs492 = DN_1_2*clhs259;
const double clhs493 = DN_3_0*clhs62;
const double clhs494 = DN_3_1*clhs194;
const double clhs495 = DN_3_2*clhs197;
const double clhs496 = DN_1_1*DN_3_1;
const double clhs497 = -clhs3*clhs496 + clhs429;
const double clhs498 = -DN_3_2*clhs451;
const double clhs499 = DN_1_0*clhs155;
const double clhs500 = DN_1_1*clhs268;
const double clhs501 = DN_1_2*clhs271;
const double clhs502 = DN_3_0*clhs66;
const double clhs503 = DN_3_1*clhs197;
const double clhs504 = DN_3_2*clhs295;
const double clhs505 = DN_1_1*N_3;
const double clhs506 = pow(DN_1_2, 2);
const double clhs507 = DN_1_0*clhs78;
const double clhs508 = DN_1_1*clhs209;
const double clhs509 = DN_1_2*clhs302;
const double clhs510 = DN_1_2*N_1;
const double clhs511 = DN_1_2*clhs3;
const double clhs512 = -DN_2_0*clhs511;
const double clhs513 = DN_1_0*clhs90;
const double clhs514 = DN_1_1*clhs219;
const double clhs515 = DN_1_2*clhs313;
const double clhs516 = DN_2_0*clhs72;
const double clhs517 = DN_2_1*clhs75;
const double clhs518 = DN_2_2*clhs78;
const double clhs519 = -DN_2_1*clhs511;
const double clhs520 = DN_1_0*clhs106;
const double clhs521 = DN_1_1*clhs228;
const double clhs522 = DN_1_2*clhs321;
const double clhs523 = DN_2_0*clhs75;
const double clhs524 = DN_2_1*clhs206;
const double clhs525 = DN_2_2*clhs209;
const double clhs526 = DN_1_0*clhs118;
const double clhs527 = DN_1_1*clhs240;
const double clhs528 = DN_1_2*clhs328;
const double clhs529 = DN_2_0*clhs78;
const double clhs530 = DN_2_1*clhs209;
const double clhs531 = DN_2_2*clhs302;
const double clhs532 = DN_1_2*DN_2_2;
const double clhs533 = -clhs3*clhs532 + clhs405;
const double clhs534 = DN_1_2*N_2;
const double clhs535 = -DN_3_0*clhs511;
const double clhs536 = DN_1_0*clhs130;
const double clhs537 = DN_1_1*clhs250;
const double clhs538 = DN_1_2*clhs339;
const double clhs539 = DN_3_0*clhs72;
const double clhs540 = DN_3_1*clhs75;
const double clhs541 = DN_3_2*clhs78;
const double clhs542 = -DN_3_1*clhs511;
const double clhs543 = DN_1_0*clhs146;
const double clhs544 = DN_1_1*clhs259;
const double clhs545 = DN_1_2*clhs347;
const double clhs546 = DN_3_0*clhs75;
const double clhs547 = DN_3_1*clhs206;
const double clhs548 = DN_3_2*clhs209;
const double clhs549 = DN_1_0*clhs158;
const double clhs550 = DN_1_1*clhs271;
const double clhs551 = DN_1_2*clhs354;
const double clhs552 = DN_3_0*clhs78;
const double clhs553 = DN_3_1*clhs209;
const double clhs554 = DN_3_2*clhs302;
const double clhs555 = DN_1_2*DN_3_2;
const double clhs556 = -clhs3*clhs555 + clhs429;
const double clhs557 = DN_1_2*N_3;
const double clhs558 = DN_2_0*N_1;
const double clhs559 = DN_2_1*N_1;
const double clhs560 = DN_2_2*N_1;
const double clhs561 = clhs364*(clhs406 + clhs473 + clhs532);
const double clhs562 = DN_3_0*N_1;
const double clhs563 = DN_3_1*N_1;
const double clhs564 = DN_3_2*N_1;
const double clhs565 = clhs364*(clhs430 + clhs496 + clhs555);
const double clhs566 = pow(N_2, 2)*clhs0;
const double clhs567 = pow(DN_2_0, 2);
const double clhs568 = DN_2_0*clhs84;
const double clhs569 = DN_2_1*clhs87;
const double clhs570 = DN_2_2*clhs90;
const double clhs571 = DN_2_0*clhs3;
const double clhs572 = -DN_2_1*clhs571;
const double clhs573 = DN_2_0*clhs99;
const double clhs574 = DN_2_1*clhs102;
const double clhs575 = DN_2_2*clhs106;
const double clhs576 = DN_2_0*clhs87;
const double clhs577 = DN_2_1*clhs217;
const double clhs578 = DN_2_2*clhs219;
const double clhs579 = -DN_2_2*clhs571;
const double clhs580 = DN_2_0*clhs112;
const double clhs581 = DN_2_1*clhs115;
const double clhs582 = DN_2_2*clhs118;
const double clhs583 = DN_2_0*clhs90;
const double clhs584 = DN_2_1*clhs219;
const double clhs585 = DN_2_2*clhs313;
const double clhs586 = DN_2_0*N_2;
const double clhs587 = DN_2_0*clhs124;
const double clhs588 = DN_2_1*clhs127;
const double clhs589 = DN_2_2*clhs130;
const double clhs590 = DN_3_0*clhs84;
const double clhs591 = DN_3_1*clhs87;
const double clhs592 = DN_3_2*clhs90;
const double clhs593 = N_2*N_3*clhs0;
const double clhs594 = DN_2_0*DN_3_0;
const double clhs595 = -clhs3*clhs594 + clhs593;
const double clhs596 = -DN_3_1*clhs571;
const double clhs597 = DN_2_0*clhs139;
const double clhs598 = DN_2_1*clhs142;
const double clhs599 = DN_2_2*clhs146;
const double clhs600 = DN_3_0*clhs87;
const double clhs601 = DN_3_1*clhs217;
const double clhs602 = DN_3_2*clhs219;
const double clhs603 = -DN_3_2*clhs571;
const double clhs604 = DN_2_0*clhs152;
const double clhs605 = DN_2_1*clhs155;
const double clhs606 = DN_2_2*clhs158;
const double clhs607 = DN_3_0*clhs90;
const double clhs608 = DN_3_1*clhs219;
const double clhs609 = DN_3_2*clhs313;
const double clhs610 = DN_2_0*N_3;
const double clhs611 = pow(DN_2_1, 2);
const double clhs612 = DN_2_0*clhs102;
const double clhs613 = DN_2_1*clhs225;
const double clhs614 = DN_2_2*clhs228;
const double clhs615 = DN_2_1*clhs3;
const double clhs616 = -DN_2_2*clhs615;
const double clhs617 = DN_2_0*clhs115;
const double clhs618 = DN_2_1*clhs237;
const double clhs619 = DN_2_2*clhs240;
const double clhs620 = DN_2_0*clhs106;
const double clhs621 = DN_2_1*clhs228;
const double clhs622 = DN_2_2*clhs321;
const double clhs623 = DN_2_1*N_2;
const double clhs624 = -DN_3_0*clhs615;
const double clhs625 = DN_2_0*clhs127;
const double clhs626 = DN_2_1*clhs248;
const double clhs627 = DN_2_2*clhs250;
const double clhs628 = DN_3_0*clhs99;
const double clhs629 = DN_3_1*clhs102;
const double clhs630 = DN_3_2*clhs106;
const double clhs631 = DN_2_0*clhs142;
const double clhs632 = DN_2_1*clhs256;
const double clhs633 = DN_2_2*clhs259;
const double clhs634 = DN_3_0*clhs102;
const double clhs635 = DN_3_1*clhs225;
const double clhs636 = DN_3_2*clhs228;
const double clhs637 = DN_2_1*DN_3_1;
const double clhs638 = -clhs3*clhs637 + clhs593;
const double clhs639 = -DN_3_2*clhs615;
const double clhs640 = DN_2_0*clhs155;
const double clhs641 = DN_2_1*clhs268;
const double clhs642 = DN_2_2*clhs271;
const double clhs643 = DN_3_0*clhs106;
const double clhs644 = DN_3_1*clhs228;
const double clhs645 = DN_3_2*clhs321;
const double clhs646 = DN_2_1*N_3;
const double clhs647 = pow(DN_2_2, 2);
const double clhs648 = DN_2_0*clhs118;
const double clhs649 = DN_2_1*clhs240;
const double clhs650 = DN_2_2*clhs328;
const double clhs651 = DN_2_2*N_2;
const double clhs652 = DN_2_2*clhs3;
const double clhs653 = -DN_3_0*clhs652;
const double clhs654 = DN_2_0*clhs130;
const double clhs655 = DN_2_1*clhs250;
const double clhs656 = DN_2_2*clhs339;
const double clhs657 = DN_3_0*clhs112;
const double clhs658 = DN_3_1*clhs115;
const double clhs659 = DN_3_2*clhs118;
const double clhs660 = -DN_3_1*clhs652;
const double clhs661 = DN_2_0*clhs146;
const double clhs662 = DN_2_1*clhs259;
const double clhs663 = DN_2_2*clhs347;
const double clhs664 = DN_3_0*clhs115;
const double clhs665 = DN_3_1*clhs237;
const double clhs666 = DN_3_2*clhs240;
const double clhs667 = DN_2_0*clhs158;
const double clhs668 = DN_2_1*clhs271;
const double clhs669 = DN_2_2*clhs354;
const double clhs670 = DN_3_0*clhs118;
const double clhs671 = DN_3_1*clhs240;
const double clhs672 = DN_3_2*clhs328;
const double clhs673 = DN_2_2*DN_3_2;
const double clhs674 = -clhs3*clhs673 + clhs593;
const double clhs675 = DN_2_2*N_3;
const double clhs676 = DN_3_0*N_2;
const double clhs677 = DN_3_1*N_2;
const double clhs678 = DN_3_2*N_2;
const double clhs679 = clhs364*(clhs594 + clhs637 + clhs673);
const double clhs680 = pow(N_3, 2)*clhs0;
const double clhs681 = pow(DN_3_0, 2);
const double clhs682 = DN_3_0*clhs124;
const double clhs683 = DN_3_1*clhs127;
const double clhs684 = DN_3_2*clhs130;
const double clhs685 = DN_3_0*clhs3;
const double clhs686 = -DN_3_1*clhs685;
const double clhs687 = DN_3_0*clhs139;
const double clhs688 = DN_3_1*clhs142;
const double clhs689 = DN_3_2*clhs146;
const double clhs690 = DN_3_0*clhs127;
const double clhs691 = DN_3_1*clhs248;
const double clhs692 = DN_3_2*clhs250;
const double clhs693 = -DN_3_2*clhs685;
const double clhs694 = DN_3_0*clhs152;
const double clhs695 = DN_3_1*clhs155;
const double clhs696 = DN_3_2*clhs158;
const double clhs697 = DN_3_0*clhs130;
const double clhs698 = DN_3_1*clhs250;
const double clhs699 = DN_3_2*clhs339;
const double clhs700 = DN_3_0*N_3;
const double clhs701 = pow(DN_3_1, 2);
const double clhs702 = DN_3_0*clhs142;
const double clhs703 = DN_3_1*clhs256;
const double clhs704 = DN_3_2*clhs259;
const double clhs705 = -DN_3_1*DN_3_2*clhs3;
const double clhs706 = DN_3_0*clhs155;
const double clhs707 = DN_3_1*clhs268;
const double clhs708 = DN_3_2*clhs271;
const double clhs709 = DN_3_0*clhs146;
const double clhs710 = DN_3_1*clhs259;
const double clhs711 = DN_3_2*clhs347;
const double clhs712 = DN_3_1*N_3;
const double clhs713 = pow(DN_3_2, 2);
const double clhs714 = DN_3_0*clhs158;
const double clhs715 = DN_3_1*clhs271;
const double clhs716 = DN_3_2*clhs354;
const double clhs717 = DN_3_2*N_3;
lhs(0,0)=clhs1 - clhs11*clhs12 + clhs11 - clhs12*clhs5 - clhs12*clhs8 - clhs2*clhs3 + clhs5 + clhs8;
lhs(0,1)=-clhs12*clhs24 - clhs12*clhs26 - clhs12*clhs29 + clhs14 + clhs16 + clhs19 + clhs23;
lhs(0,2)=-clhs12*clhs38 - clhs12*clhs39 - clhs12*clhs41 + clhs30 + clhs32 + clhs34 + clhs37;
lhs(0,3)=-clhs42;
lhs(0,4)=-clhs12*clhs51 - clhs12*clhs52 - clhs12*clhs53 + clhs44 + clhs47 + clhs50 + clhs57;
lhs(0,5)=-clhs12*clhs68 - clhs12*clhs69 - clhs12*clhs70 + clhs58 + clhs60 + clhs63 + clhs67;
lhs(0,6)=-clhs12*clhs80 - clhs12*clhs81 - clhs12*clhs82 + clhs71 + clhs73 + clhs76 + clhs79;
lhs(0,7)=-clhs83;
lhs(0,8)=-clhs12*clhs92 - clhs12*clhs93 - clhs12*clhs94 + clhs85 + clhs88 + clhs91 + clhs97;
lhs(0,9)=clhs100 + clhs103 + clhs107 - clhs108*clhs12 - clhs109*clhs12 - clhs110*clhs12 + clhs98;
lhs(0,10)=clhs111 + clhs113 + clhs116 + clhs119 - clhs12*clhs120 - clhs12*clhs121 - clhs12*clhs122;
lhs(0,11)=-clhs123;
lhs(0,12)=-clhs12*clhs132 - clhs12*clhs133 - clhs12*clhs134 + clhs125 + clhs128 + clhs131 + clhs137;
lhs(0,13)=-clhs12*clhs148 - clhs12*clhs149 - clhs12*clhs150 + clhs138 + clhs140 + clhs143 + clhs147;
lhs(0,14)=-clhs12*clhs160 - clhs12*clhs161 - clhs12*clhs162 + clhs151 + clhs153 + clhs156 + clhs159;
lhs(0,15)=-clhs163;
lhs(1,0)=-clhs12*clhs16 - clhs12*clhs19 - clhs12*clhs23 + clhs14 + clhs24 + clhs26 + clhs29;
lhs(1,1)=clhs1 - clhs12*clhs165 - clhs12*clhs167 - clhs12*clhs170 - clhs164*clhs3 + clhs165 + clhs167 + clhs170;
lhs(1,2)=-clhs12*clhs179 - clhs12*clhs180 - clhs12*clhs182 + clhs172 + clhs173 + clhs175 + clhs178;
lhs(1,3)=-clhs183;
lhs(1,4)=-clhs12*clhs190 - clhs12*clhs191 - clhs12*clhs192 + clhs184 + clhs185 + clhs187 + clhs189;
lhs(1,5)=-clhs12*clhs199 - clhs12*clhs200 - clhs12*clhs201 + clhs193 + clhs195 + clhs198 + clhs203;
lhs(1,6)=-clhs12*clhs211 - clhs12*clhs212 - clhs12*clhs213 + clhs204 + clhs205 + clhs207 + clhs210;
lhs(1,7)=-clhs214;
lhs(1,8)=-clhs12*clhs221 - clhs12*clhs222 - clhs12*clhs223 + clhs215 + clhs216 + clhs218 + clhs220;
lhs(1,9)=-clhs12*clhs230 - clhs12*clhs231 - clhs12*clhs232 + clhs224 + clhs226 + clhs229 + clhs234;
lhs(1,10)=-clhs12*clhs242 - clhs12*clhs243 - clhs12*clhs244 + clhs235 + clhs236 + clhs238 + clhs241;
lhs(1,11)=-clhs245;
lhs(1,12)=-clhs12*clhs252 - clhs12*clhs253 - clhs12*clhs254 + clhs246 + clhs247 + clhs249 + clhs251;
lhs(1,13)=-clhs12*clhs261 - clhs12*clhs262 - clhs12*clhs263 + clhs255 + clhs257 + clhs260 + clhs265;
lhs(1,14)=-clhs12*clhs273 - clhs12*clhs274 - clhs12*clhs275 + clhs266 + clhs267 + clhs269 + clhs272;
lhs(1,15)=-clhs276;
lhs(2,0)=-clhs12*clhs32 - clhs12*clhs34 - clhs12*clhs37 + clhs30 + clhs38 + clhs39 + clhs41;
lhs(2,1)=-clhs12*clhs173 - clhs12*clhs175 - clhs12*clhs178 + clhs172 + clhs179 + clhs180 + clhs182;
lhs(2,2)=clhs1 - clhs12*clhs278 - clhs12*clhs279 - clhs12*clhs281 - clhs277*clhs3 + clhs278 + clhs279 + clhs281;
lhs(2,3)=-clhs282;
lhs(2,4)=-clhs12*clhs289 - clhs12*clhs290 - clhs12*clhs291 + clhs284 + clhs285 + clhs286 + clhs288;
lhs(2,5)=-clhs12*clhs297 - clhs12*clhs298 - clhs12*clhs299 + clhs292 + clhs293 + clhs294 + clhs296;
lhs(2,6)=-clhs12*clhs304 - clhs12*clhs305 - clhs12*clhs306 + clhs300 + clhs301 + clhs303 + clhs308;
lhs(2,7)=-clhs309;
lhs(2,8)=-clhs12*clhs315 - clhs12*clhs316 - clhs12*clhs317 + clhs310 + clhs311 + clhs312 + clhs314;
lhs(2,9)=-clhs12*clhs323 - clhs12*clhs324 - clhs12*clhs325 + clhs318 + clhs319 + clhs320 + clhs322;
lhs(2,10)=-clhs12*clhs330 - clhs12*clhs331 - clhs12*clhs332 + clhs326 + clhs327 + clhs329 + clhs334;
lhs(2,11)=-clhs335;
lhs(2,12)=-clhs12*clhs341 - clhs12*clhs342 - clhs12*clhs343 + clhs336 + clhs337 + clhs338 + clhs340;
lhs(2,13)=-clhs12*clhs349 - clhs12*clhs350 - clhs12*clhs351 + clhs344 + clhs345 + clhs346 + clhs348;
lhs(2,14)=-clhs12*clhs356 - clhs12*clhs357 - clhs12*clhs358 + clhs352 + clhs353 + clhs355 + clhs360;
lhs(2,15)=-clhs361;
lhs(3,0)=clhs363*clhs42;
lhs(3,1)=clhs183*clhs363;
lhs(3,2)=clhs282*clhs363;
lhs(3,3)=clhs364*(clhs164 + clhs2 + clhs277);
lhs(3,4)=clhs362*clhs83 + clhs365;
lhs(3,5)=clhs214*clhs362 + clhs366;
lhs(3,6)=clhs309*clhs362 + clhs367;
lhs(3,7)=clhs368;
lhs(3,8)=clhs123*clhs362 + clhs369;
lhs(3,9)=clhs245*clhs362 + clhs370;
lhs(3,10)=clhs335*clhs362 + clhs371;
lhs(3,11)=clhs372;
lhs(3,12)=clhs163*clhs362 + clhs373;
lhs(3,13)=clhs276*clhs362 + clhs374;
lhs(3,14)=clhs361*clhs362 + clhs375;
lhs(3,15)=clhs376;
lhs(4,0)=-clhs12*clhs44 - clhs12*clhs47 - clhs12*clhs50 + clhs51 + clhs52 + clhs53 + clhs57;
lhs(4,1)=-clhs12*clhs185 - clhs12*clhs187 - clhs12*clhs189 + clhs184 + clhs190 + clhs191 + clhs192;
lhs(4,2)=-clhs12*clhs285 - clhs12*clhs286 - clhs12*clhs288 + clhs284 + clhs289 + clhs290 + clhs291;
lhs(4,3)=-clhs365;
lhs(4,4)=-clhs12*clhs379 - clhs12*clhs380 - clhs12*clhs381 - clhs3*clhs378 + clhs377 + clhs379 + clhs380 + clhs381;
lhs(4,5)=-clhs12*clhs387 - clhs12*clhs388 - clhs12*clhs389 + clhs383 + clhs384 + clhs385 + clhs386;
lhs(4,6)=-clhs12*clhs394 - clhs12*clhs395 - clhs12*clhs396 + clhs390 + clhs391 + clhs392 + clhs393;
lhs(4,7)=-clhs397;
lhs(4,8)=-clhs12*clhs401 - clhs12*clhs402 - clhs12*clhs403 + clhs398 + clhs399 + clhs400 + clhs407;
lhs(4,9)=-clhs12*clhs412 - clhs12*clhs413 - clhs12*clhs414 + clhs408 + clhs409 + clhs410 + clhs411;
lhs(4,10)=-clhs12*clhs419 - clhs12*clhs420 - clhs12*clhs421 + clhs415 + clhs416 + clhs417 + clhs418;
lhs(4,11)=-clhs422;
lhs(4,12)=-clhs12*clhs426 - clhs12*clhs427 - clhs12*clhs428 + clhs423 + clhs424 + clhs425 + clhs431;
lhs(4,13)=-clhs12*clhs436 - clhs12*clhs437 - clhs12*clhs438 + clhs432 + clhs433 + clhs434 + clhs435;
lhs(4,14)=-clhs12*clhs443 - clhs12*clhs444 - clhs12*clhs445 + clhs439 + clhs440 + clhs441 + clhs442;
lhs(4,15)=-clhs446;
lhs(5,0)=-clhs12*clhs60 - clhs12*clhs63 - clhs12*clhs67 + clhs58 + clhs68 + clhs69 + clhs70;
lhs(5,1)=-clhs12*clhs193 - clhs12*clhs195 - clhs12*clhs198 + clhs199 + clhs200 + clhs201 + clhs203;
lhs(5,2)=-clhs12*clhs293 - clhs12*clhs294 - clhs12*clhs296 + clhs292 + clhs297 + clhs298 + clhs299;
lhs(5,3)=-clhs366;
lhs(5,4)=-clhs12*clhs384 - clhs12*clhs385 - clhs12*clhs386 + clhs383 + clhs387 + clhs388 + clhs389;
lhs(5,5)=-clhs12*clhs448 - clhs12*clhs449 - clhs12*clhs450 - clhs3*clhs447 + clhs377 + clhs448 + clhs449 + clhs450;
lhs(5,6)=-clhs12*clhs456 - clhs12*clhs457 - clhs12*clhs458 + clhs452 + clhs453 + clhs454 + clhs455;
lhs(5,7)=-clhs459;
lhs(5,8)=-clhs12*clhs464 - clhs12*clhs465 - clhs12*clhs466 + clhs460 + clhs461 + clhs462 + clhs463;
lhs(5,9)=-clhs12*clhs470 - clhs12*clhs471 - clhs12*clhs472 + clhs467 + clhs468 + clhs469 + clhs474;
lhs(5,10)=-clhs12*clhs479 - clhs12*clhs480 - clhs12*clhs481 + clhs475 + clhs476 + clhs477 + clhs478;
lhs(5,11)=-clhs482;
lhs(5,12)=-clhs12*clhs487 - clhs12*clhs488 - clhs12*clhs489 + clhs483 + clhs484 + clhs485 + clhs486;
lhs(5,13)=-clhs12*clhs493 - clhs12*clhs494 - clhs12*clhs495 + clhs490 + clhs491 + clhs492 + clhs497;
lhs(5,14)=-clhs12*clhs502 - clhs12*clhs503 - clhs12*clhs504 + clhs498 + clhs499 + clhs500 + clhs501;
lhs(5,15)=-clhs505;
lhs(6,0)=-clhs12*clhs73 - clhs12*clhs76 - clhs12*clhs79 + clhs71 + clhs80 + clhs81 + clhs82;
lhs(6,1)=-clhs12*clhs205 - clhs12*clhs207 - clhs12*clhs210 + clhs204 + clhs211 + clhs212 + clhs213;
lhs(6,2)=-clhs12*clhs300 - clhs12*clhs301 - clhs12*clhs303 + clhs304 + clhs305 + clhs306 + clhs308;
lhs(6,3)=-clhs367;
lhs(6,4)=-clhs12*clhs391 - clhs12*clhs392 - clhs12*clhs393 + clhs390 + clhs394 + clhs395 + clhs396;
lhs(6,5)=-clhs12*clhs453 - clhs12*clhs454 - clhs12*clhs455 + clhs452 + clhs456 + clhs457 + clhs458;
lhs(6,6)=-clhs12*clhs507 - clhs12*clhs508 - clhs12*clhs509 - clhs3*clhs506 + clhs377 + clhs507 + clhs508 + clhs509;
lhs(6,7)=-clhs510;
lhs(6,8)=-clhs12*clhs516 - clhs12*clhs517 - clhs12*clhs518 + clhs512 + clhs513 + clhs514 + clhs515;
lhs(6,9)=-clhs12*clhs523 - clhs12*clhs524 - clhs12*clhs525 + clhs519 + clhs520 + clhs521 + clhs522;
lhs(6,10)=-clhs12*clhs529 - clhs12*clhs530 - clhs12*clhs531 + clhs526 + clhs527 + clhs528 + clhs533;
lhs(6,11)=-clhs534;
lhs(6,12)=-clhs12*clhs539 - clhs12*clhs540 - clhs12*clhs541 + clhs535 + clhs536 + clhs537 + clhs538;
lhs(6,13)=-clhs12*clhs546 - clhs12*clhs547 - clhs12*clhs548 + clhs542 + clhs543 + clhs544 + clhs545;
lhs(6,14)=-clhs12*clhs552 - clhs12*clhs553 - clhs12*clhs554 + clhs549 + clhs550 + clhs551 + clhs556;
lhs(6,15)=-clhs557;
lhs(7,0)=clhs362*clhs365 + clhs83;
lhs(7,1)=clhs214 + clhs362*clhs366;
lhs(7,2)=clhs309 + clhs362*clhs367;
lhs(7,3)=clhs368;
lhs(7,4)=clhs363*clhs397;
lhs(7,5)=clhs363*clhs459;
lhs(7,6)=clhs363*clhs510;
lhs(7,7)=clhs364*(clhs378 + clhs447 + clhs506);
lhs(7,8)=clhs362*clhs422 + clhs558;
lhs(7,9)=clhs362*clhs482 + clhs559;
lhs(7,10)=clhs362*clhs534 + clhs560;
lhs(7,11)=clhs561;
lhs(7,12)=clhs362*clhs446 + clhs562;
lhs(7,13)=clhs362*clhs505 + clhs563;
lhs(7,14)=clhs362*clhs557 + clhs564;
lhs(7,15)=clhs565;
lhs(8,0)=-clhs12*clhs85 - clhs12*clhs88 - clhs12*clhs91 + clhs92 + clhs93 + clhs94 + clhs97;
lhs(8,1)=-clhs12*clhs216 - clhs12*clhs218 - clhs12*clhs220 + clhs215 + clhs221 + clhs222 + clhs223;
lhs(8,2)=-clhs12*clhs311 - clhs12*clhs312 - clhs12*clhs314 + clhs310 + clhs315 + clhs316 + clhs317;
lhs(8,3)=-clhs369;
lhs(8,4)=-clhs12*clhs398 - clhs12*clhs399 - clhs12*clhs400 + clhs401 + clhs402 + clhs403 + clhs407;
lhs(8,5)=-clhs12*clhs461 - clhs12*clhs462 - clhs12*clhs463 + clhs460 + clhs464 + clhs465 + clhs466;
lhs(8,6)=-clhs12*clhs513 - clhs12*clhs514 - clhs12*clhs515 + clhs512 + clhs516 + clhs517 + clhs518;
lhs(8,7)=-clhs558;
lhs(8,8)=-clhs12*clhs568 - clhs12*clhs569 - clhs12*clhs570 - clhs3*clhs567 + clhs566 + clhs568 + clhs569 + clhs570;
lhs(8,9)=-clhs12*clhs576 - clhs12*clhs577 - clhs12*clhs578 + clhs572 + clhs573 + clhs574 + clhs575;
lhs(8,10)=-clhs12*clhs583 - clhs12*clhs584 - clhs12*clhs585 + clhs579 + clhs580 + clhs581 + clhs582;
lhs(8,11)=-clhs586;
lhs(8,12)=-clhs12*clhs590 - clhs12*clhs591 - clhs12*clhs592 + clhs587 + clhs588 + clhs589 + clhs595;
lhs(8,13)=-clhs12*clhs600 - clhs12*clhs601 - clhs12*clhs602 + clhs596 + clhs597 + clhs598 + clhs599;
lhs(8,14)=-clhs12*clhs607 - clhs12*clhs608 - clhs12*clhs609 + clhs603 + clhs604 + clhs605 + clhs606;
lhs(8,15)=-clhs610;
lhs(9,0)=-clhs100*clhs12 - clhs103*clhs12 - clhs107*clhs12 + clhs108 + clhs109 + clhs110 + clhs98;
lhs(9,1)=-clhs12*clhs224 - clhs12*clhs226 - clhs12*clhs229 + clhs230 + clhs231 + clhs232 + clhs234;
lhs(9,2)=-clhs12*clhs319 - clhs12*clhs320 - clhs12*clhs322 + clhs318 + clhs323 + clhs324 + clhs325;
lhs(9,3)=-clhs370;
lhs(9,4)=-clhs12*clhs409 - clhs12*clhs410 - clhs12*clhs411 + clhs408 + clhs412 + clhs413 + clhs414;
lhs(9,5)=-clhs12*clhs467 - clhs12*clhs468 - clhs12*clhs469 + clhs470 + clhs471 + clhs472 + clhs474;
lhs(9,6)=-clhs12*clhs520 - clhs12*clhs521 - clhs12*clhs522 + clhs519 + clhs523 + clhs524 + clhs525;
lhs(9,7)=-clhs559;
lhs(9,8)=-clhs12*clhs573 - clhs12*clhs574 - clhs12*clhs575 + clhs572 + clhs576 + clhs577 + clhs578;
lhs(9,9)=-clhs12*clhs612 - clhs12*clhs613 - clhs12*clhs614 - clhs3*clhs611 + clhs566 + clhs612 + clhs613 + clhs614;
lhs(9,10)=-clhs12*clhs620 - clhs12*clhs621 - clhs12*clhs622 + clhs616 + clhs617 + clhs618 + clhs619;
lhs(9,11)=-clhs623;
lhs(9,12)=-clhs12*clhs628 - clhs12*clhs629 - clhs12*clhs630 + clhs624 + clhs625 + clhs626 + clhs627;
lhs(9,13)=-clhs12*clhs634 - clhs12*clhs635 - clhs12*clhs636 + clhs631 + clhs632 + clhs633 + clhs638;
lhs(9,14)=-clhs12*clhs643 - clhs12*clhs644 - clhs12*clhs645 + clhs639 + clhs640 + clhs641 + clhs642;
lhs(9,15)=-clhs646;
lhs(10,0)=clhs111 - clhs113*clhs12 - clhs116*clhs12 - clhs119*clhs12 + clhs120 + clhs121 + clhs122;
lhs(10,1)=-clhs12*clhs236 - clhs12*clhs238 - clhs12*clhs241 + clhs235 + clhs242 + clhs243 + clhs244;
lhs(10,2)=-clhs12*clhs326 - clhs12*clhs327 - clhs12*clhs329 + clhs330 + clhs331 + clhs332 + clhs334;
lhs(10,3)=-clhs371;
lhs(10,4)=-clhs12*clhs416 - clhs12*clhs417 - clhs12*clhs418 + clhs415 + clhs419 + clhs420 + clhs421;
lhs(10,5)=-clhs12*clhs476 - clhs12*clhs477 - clhs12*clhs478 + clhs475 + clhs479 + clhs480 + clhs481;
lhs(10,6)=-clhs12*clhs526 - clhs12*clhs527 - clhs12*clhs528 + clhs529 + clhs530 + clhs531 + clhs533;
lhs(10,7)=-clhs560;
lhs(10,8)=-clhs12*clhs580 - clhs12*clhs581 - clhs12*clhs582 + clhs579 + clhs583 + clhs584 + clhs585;
lhs(10,9)=-clhs12*clhs617 - clhs12*clhs618 - clhs12*clhs619 + clhs616 + clhs620 + clhs621 + clhs622;
lhs(10,10)=-clhs12*clhs648 - clhs12*clhs649 - clhs12*clhs650 - clhs3*clhs647 + clhs566 + clhs648 + clhs649 + clhs650;
lhs(10,11)=-clhs651;
lhs(10,12)=-clhs12*clhs657 - clhs12*clhs658 - clhs12*clhs659 + clhs653 + clhs654 + clhs655 + clhs656;
lhs(10,13)=-clhs12*clhs664 - clhs12*clhs665 - clhs12*clhs666 + clhs660 + clhs661 + clhs662 + clhs663;
lhs(10,14)=-clhs12*clhs670 - clhs12*clhs671 - clhs12*clhs672 + clhs667 + clhs668 + clhs669 + clhs674;
lhs(10,15)=-clhs675;
lhs(11,0)=clhs123 + clhs362*clhs369;
lhs(11,1)=clhs245 + clhs362*clhs370;
lhs(11,2)=clhs335 + clhs362*clhs371;
lhs(11,3)=clhs372;
lhs(11,4)=clhs362*clhs558 + clhs422;
lhs(11,5)=clhs362*clhs559 + clhs482;
lhs(11,6)=clhs362*clhs560 + clhs534;
lhs(11,7)=clhs561;
lhs(11,8)=clhs363*clhs586;
lhs(11,9)=clhs363*clhs623;
lhs(11,10)=clhs363*clhs651;
lhs(11,11)=clhs364*(clhs567 + clhs611 + clhs647);
lhs(11,12)=clhs362*clhs610 + clhs676;
lhs(11,13)=clhs362*clhs646 + clhs677;
lhs(11,14)=clhs362*clhs675 + clhs678;
lhs(11,15)=clhs679;
lhs(12,0)=-clhs12*clhs125 - clhs12*clhs128 - clhs12*clhs131 + clhs132 + clhs133 + clhs134 + clhs137;
lhs(12,1)=-clhs12*clhs247 - clhs12*clhs249 - clhs12*clhs251 + clhs246 + clhs252 + clhs253 + clhs254;
lhs(12,2)=-clhs12*clhs337 - clhs12*clhs338 - clhs12*clhs340 + clhs336 + clhs341 + clhs342 + clhs343;
lhs(12,3)=-clhs373;
lhs(12,4)=-clhs12*clhs423 - clhs12*clhs424 - clhs12*clhs425 + clhs426 + clhs427 + clhs428 + clhs431;
lhs(12,5)=-clhs12*clhs484 - clhs12*clhs485 - clhs12*clhs486 + clhs483 + clhs487 + clhs488 + clhs489;
lhs(12,6)=-clhs12*clhs536 - clhs12*clhs537 - clhs12*clhs538 + clhs535 + clhs539 + clhs540 + clhs541;
lhs(12,7)=-clhs562;
lhs(12,8)=-clhs12*clhs587 - clhs12*clhs588 - clhs12*clhs589 + clhs590 + clhs591 + clhs592 + clhs595;
lhs(12,9)=-clhs12*clhs625 - clhs12*clhs626 - clhs12*clhs627 + clhs624 + clhs628 + clhs629 + clhs630;
lhs(12,10)=-clhs12*clhs654 - clhs12*clhs655 - clhs12*clhs656 + clhs653 + clhs657 + clhs658 + clhs659;
lhs(12,11)=-clhs676;
lhs(12,12)=-clhs12*clhs682 - clhs12*clhs683 - clhs12*clhs684 - clhs3*clhs681 + clhs680 + clhs682 + clhs683 + clhs684;
lhs(12,13)=-clhs12*clhs690 - clhs12*clhs691 - clhs12*clhs692 + clhs686 + clhs687 + clhs688 + clhs689;
lhs(12,14)=-clhs12*clhs697 - clhs12*clhs698 - clhs12*clhs699 + clhs693 + clhs694 + clhs695 + clhs696;
lhs(12,15)=-clhs700;
lhs(13,0)=-clhs12*clhs140 - clhs12*clhs143 - clhs12*clhs147 + clhs138 + clhs148 + clhs149 + clhs150;
lhs(13,1)=-clhs12*clhs255 - clhs12*clhs257 - clhs12*clhs260 + clhs261 + clhs262 + clhs263 + clhs265;
lhs(13,2)=-clhs12*clhs345 - clhs12*clhs346 - clhs12*clhs348 + clhs344 + clhs349 + clhs350 + clhs351;
lhs(13,3)=-clhs374;
lhs(13,4)=-clhs12*clhs433 - clhs12*clhs434 - clhs12*clhs435 + clhs432 + clhs436 + clhs437 + clhs438;
lhs(13,5)=-clhs12*clhs490 - clhs12*clhs491 - clhs12*clhs492 + clhs493 + clhs494 + clhs495 + clhs497;
lhs(13,6)=-clhs12*clhs543 - clhs12*clhs544 - clhs12*clhs545 + clhs542 + clhs546 + clhs547 + clhs548;
lhs(13,7)=-clhs563;
lhs(13,8)=-clhs12*clhs597 - clhs12*clhs598 - clhs12*clhs599 + clhs596 + clhs600 + clhs601 + clhs602;
lhs(13,9)=-clhs12*clhs631 - clhs12*clhs632 - clhs12*clhs633 + clhs634 + clhs635 + clhs636 + clhs638;
lhs(13,10)=-clhs12*clhs661 - clhs12*clhs662 - clhs12*clhs663 + clhs660 + clhs664 + clhs665 + clhs666;
lhs(13,11)=-clhs677;
lhs(13,12)=-clhs12*clhs687 - clhs12*clhs688 - clhs12*clhs689 + clhs686 + clhs690 + clhs691 + clhs692;
lhs(13,13)=-clhs12*clhs702 - clhs12*clhs703 - clhs12*clhs704 - clhs3*clhs701 + clhs680 + clhs702 + clhs703 + clhs704;
lhs(13,14)=-clhs12*clhs709 - clhs12*clhs710 - clhs12*clhs711 + clhs705 + clhs706 + clhs707 + clhs708;
lhs(13,15)=-clhs712;
lhs(14,0)=-clhs12*clhs153 - clhs12*clhs156 - clhs12*clhs159 + clhs151 + clhs160 + clhs161 + clhs162;
lhs(14,1)=-clhs12*clhs267 - clhs12*clhs269 - clhs12*clhs272 + clhs266 + clhs273 + clhs274 + clhs275;
lhs(14,2)=-clhs12*clhs352 - clhs12*clhs353 - clhs12*clhs355 + clhs356 + clhs357 + clhs358 + clhs360;
lhs(14,3)=-clhs375;
lhs(14,4)=-clhs12*clhs440 - clhs12*clhs441 - clhs12*clhs442 + clhs439 + clhs443 + clhs444 + clhs445;
lhs(14,5)=-clhs12*clhs499 - clhs12*clhs500 - clhs12*clhs501 + clhs498 + clhs502 + clhs503 + clhs504;
lhs(14,6)=-clhs12*clhs549 - clhs12*clhs550 - clhs12*clhs551 + clhs552 + clhs553 + clhs554 + clhs556;
lhs(14,7)=-clhs564;
lhs(14,8)=-clhs12*clhs604 - clhs12*clhs605 - clhs12*clhs606 + clhs603 + clhs607 + clhs608 + clhs609;
lhs(14,9)=-clhs12*clhs640 - clhs12*clhs641 - clhs12*clhs642 + clhs639 + clhs643 + clhs644 + clhs645;
lhs(14,10)=-clhs12*clhs667 - clhs12*clhs668 - clhs12*clhs669 + clhs670 + clhs671 + clhs672 + clhs674;
lhs(14,11)=-clhs678;
lhs(14,12)=-clhs12*clhs694 - clhs12*clhs695 - clhs12*clhs696 + clhs693 + clhs697 + clhs698 + clhs699;
lhs(14,13)=-clhs12*clhs706 - clhs12*clhs707 - clhs12*clhs708 + clhs705 + clhs709 + clhs710 + clhs711;
lhs(14,14)=-clhs12*clhs714 - clhs12*clhs715 - clhs12*clhs716 - clhs3*clhs713 + clhs680 + clhs714 + clhs715 + clhs716;
lhs(14,15)=-clhs717;
lhs(15,0)=clhs163 + clhs362*clhs373;
lhs(15,1)=clhs276 + clhs362*clhs374;
lhs(15,2)=clhs361 + clhs362*clhs375;
lhs(15,3)=clhs376;
lhs(15,4)=clhs362*clhs562 + clhs446;
lhs(15,5)=clhs362*clhs563 + clhs505;
lhs(15,6)=clhs362*clhs564 + clhs557;
lhs(15,7)=clhs565;
lhs(15,8)=clhs362*clhs676 + clhs610;
lhs(15,9)=clhs362*clhs677 + clhs646;
lhs(15,10)=clhs362*clhs678 + clhs675;
lhs(15,11)=clhs679;
lhs(15,12)=clhs363*clhs700;
lhs(15,13)=clhs363*clhs712;
lhs(15,14)=clhs363*clhs717;
lhs(15,15)=clhs364*(clhs681 + clhs701 + clhs713);


    }

void Stokes3DTwoFluid::ComputeGaussPointRHSContribution(array_1d<double,16>& rhs, const element_data<4,3>& data)
    {
        const int nnodes = 4;
        const int dim = 3;

        const double rho = inner_prod(data.N, data.rho);
        const double& bdf0 = data.bdf0;
        const double& bdf1 = data.bdf1;
        const double& bdf2 = data.bdf2;

        const BoundedMatrix<double,nnodes,dim>& v = data.v;
        const BoundedMatrix<double,nnodes,dim>& vn = data.vn;
        const BoundedMatrix<double,nnodes,dim>& vnn = data.vnn;
        const BoundedMatrix<double,nnodes,dim>& f = data.f;
        const array_1d<double,nnodes>& p = data.p;

        //get constitutive matrix
        const Matrix& C = data.C;
        const Vector& stress = data.stress;

        //get shape function values
        const BoundedMatrix<double,nnodes,dim>& DN = data.DN_DX;
        const array_1d<double,nnodes>& N = data.N;

        //compute an equivalent tau by Bitrans*c*Bi
        const double tau_denom = (C(3,3) + C(4,4) + C(5,5))*(pow(DN(0,0), 2) + pow(DN(0,1), 2) + pow(DN(0,2), 2) + pow(DN(1,0), 2) + pow(DN(1,1), 2) + pow(DN(1,2), 2) + pow(DN(2,0), 2) + pow(DN(2,1), 2) + pow(DN(2,2), 2) + pow(DN(3,0), 2) + pow(DN(3,1), 2) + pow(DN(3,2), 2));

        const double tau1 = 1.0/(tau_denom*rho);
        const double tau2 = (C(3,3) + C(4,4) + C(5,5))/(6.0*rho);

        //auxiliary variables used in the calculation of the RHS
        const array_1d<double,dim> fgauss = prod(trans(f), N);
        const array_1d<double,dim> vgauss = prod(trans(v), N);
        const array_1d<double,dim> grad_p = prod(trans(DN), p);
        const double pgauss = inner_prod(N,p);

        array_1d<double,dim> acch = bdf0*vgauss;
        noalias(acch) += bdf1*prod(trans(vn), N);
        noalias(acch) += bdf2*prod(trans(vnn), N);

        const double fgauss_0 = fgauss[0];
        const double fgauss_1 = fgauss[1];
        const double fgauss_2 = fgauss[2];
        
        const double grad_p_0 = grad_p[0];
        const double grad_p_1 = grad_p[1];
        const double grad_p_2 = grad_p[2];

        const double stress_0 = stress[0];
        const double stress_1 = stress[1];
        const double stress_2 = stress[2];
        const double stress_3 = stress[3];
        const double stress_4 = stress[4];
        const double stress_5 = stress[5];

        const double N_0 = N[0];
        const double N_1 = N[1];
        const double N_2 = N[2];
        const double N_3 = N[3];

        const double f_0_0 = f(0,0); const double f_0_1 = f(0,1); const double f_0_2 = f(0,2);
        const double f_1_0 = f(0,0); const double f_1_1 = f(0,1); const double f_1_2 = f(0,2);
        const double f_2_0 = f(0,0); const double f_2_1 = f(0,1); const double f_2_2 = f(0,2);
        const double f_3_0 = f(0,0); const double f_3_1 = f(0,1); const double f_3_2 = f(0,2);

        const double vn_0_0 = vn(0,0); const double vn_0_1 = vn(0,1); const double vn_0_2 = vn(0,2);
        const double vn_1_0 = vn(1,0); const double vn_1_1 = vn(1,1); const double vn_1_2 = vn(1,2);
        const double vn_2_0 = vn(2,0); const double vn_2_1 = vn(2,1); const double vn_2_2 = vn(2,2);
        const double vn_3_0 = vn(3,0); const double vn_3_1 = vn(3,1); const double vn_3_2 = vn(3,2);

        const double vnn_0_0 = vnn(0,0); const double vnn_0_1 = vnn(0,1); const double vnn_0_2 = vnn(0,2);
        const double vnn_1_0 = vnn(1,0); const double vnn_1_1 = vnn(1,1); const double vnn_1_2 = vnn(1,2);
        const double vnn_2_0 = vnn(2,0); const double vnn_2_1 = vnn(2,1); const double vnn_2_2 = vnn(2,2);
        const double vnn_3_0 = vnn(3,0); const double vnn_3_1 = vnn(3,1); const double vnn_3_2 = vnn(3,2);

        const double acch_0 = acch[0];
        const double acch_1 = acch[1];
        const double acch_2 = acch[2];

        const double C_0_0 = C(0,0); const double C_0_1 = C(0,1); const double C_0_2 = C(0,2); const double C_0_3 = C(0,3); const double C_0_4 = C(0,4); const double C_0_5 = C(0,5);
        const double C_1_1 = C(1,1); const double C_1_2 = C(1,2); const double C_1_3 = C(1,3); const double C_1_4 = C(1,4); const double C_1_5 = C(1,5);
        const double C_2_2 = C(2,2); const double C_2_3 = C(2,3); const double C_2_4 = C(2,4); const double C_2_5 = C(2,5);
        const double C_3_3 = C(3,3); const double C_3_4 = C(3,4); const double C_3_5 = C(3,5);
        const double C_4_4 = C(4,4); const double C_4_5 = C(4,5); const double C_5_5 = C(5,5);
        
        const double DN_0_0 = DN(0,0); const double DN_0_1 = DN(0,1); const double DN_0_2 = DN(0,2);
        const double DN_1_0 = DN(1,0); const double DN_1_1 = DN(1,1); const double DN_1_2 = DN(1,2);
        const double DN_2_0 = DN(2,0); const double DN_2_1 = DN(2,1); const double DN_2_2 = DN(2,2);
        const double DN_3_0 = DN(3,0); const double DN_3_1 = DN(3,1); const double DN_3_2 = DN(3,2);
        
        const double v_0_0 = v(0,0); const double v_0_1 = v(0,1); const double v_0_2 = v(0,2);
        const double v_1_0 = v(1,0); const double v_1_1 = v(1,1); const double v_1_2 = v(1,2);
        const double v_2_0 = v(2,0); const double v_2_1 = v(2,1); const double v_2_2 = v(2,2);
        const double v_3_0 = v(3,0); const double v_3_1 = v(3,1); const double v_3_2 = v(3,2);

        const double crhs0 = acch_0*rho - fgauss_0;
const double crhs1 = DN_0_0*rho;
const double crhs2 = DN_0_0*v_0_0 + DN_0_1*v_0_1 + DN_0_2*v_0_2 + DN_1_0*v_1_0 + DN_1_1*v_1_1 + DN_1_2*v_1_2 + DN_2_0*v_2_0 + DN_2_1*v_2_1 + DN_2_2*v_2_2 + DN_3_0*v_3_0 + DN_3_1*v_3_1 + DN_3_2*v_3_2;
const double crhs3 = crhs2*tau2;
const double crhs4 = f_0_0 - rho*(bdf0*v_0_0 + bdf1*vn_0_0 + bdf2*vnn_0_0);
const double crhs5 = f_1_0 - rho*(bdf0*v_1_0 + bdf1*vn_1_0 + bdf2*vnn_1_0);
const double crhs6 = f_2_0 - rho*(bdf0*v_2_0 + bdf1*vn_2_0 + bdf2*vnn_2_0);
const double crhs7 = f_3_0 - rho*(bdf0*v_3_0 + bdf1*vn_3_0 + bdf2*vnn_3_0);
const double crhs8 = tau1*(DN_0_0*crhs4 + DN_1_0*crhs5 + DN_2_0*crhs6 + DN_3_0*crhs7);
const double crhs9 = C_1_3*DN_0_1;
const double crhs10 = f_0_1 - rho*(bdf0*v_0_1 + bdf1*vn_0_1 + bdf2*vnn_0_1);
const double crhs11 = f_1_1 - rho*(bdf0*v_1_1 + bdf1*vn_1_1 + bdf2*vnn_1_1);
const double crhs12 = f_2_1 - rho*(bdf0*v_2_1 + bdf1*vn_2_1 + bdf2*vnn_2_1);
const double crhs13 = f_3_1 - rho*(bdf0*v_3_1 + bdf1*vn_3_1 + bdf2*vnn_3_1);
const double crhs14 = tau1*(DN_0_1*crhs10 + DN_1_1*crhs11 + DN_2_1*crhs12 + DN_3_1*crhs13);
const double crhs15 = C_2_5*DN_0_2;
const double crhs16 = f_0_2 - rho*(bdf0*v_0_2 + bdf1*vn_0_2 + bdf2*vnn_0_2);
const double crhs17 = f_1_2 - rho*(bdf0*v_1_2 + bdf1*vn_1_2 + bdf2*vnn_1_2);
const double crhs18 = f_2_2 - rho*(bdf0*v_2_2 + bdf1*vn_2_2 + bdf2*vnn_2_2);
const double crhs19 = f_3_2 - rho*(bdf0*v_3_2 + bdf1*vn_3_2 + bdf2*vnn_3_2);
const double crhs20 = tau1*(DN_0_2*crhs16 + DN_1_2*crhs17 + DN_2_2*crhs18 + DN_3_2*crhs19);
const double crhs21 = C_0_3*DN_0_0;
const double crhs22 = tau1*(DN_0_0*crhs10 + DN_0_1*crhs4 + DN_1_0*crhs11 + DN_1_1*crhs5 + DN_2_0*crhs12 + DN_2_1*crhs6 + DN_3_0*crhs13 + DN_3_1*crhs7);
const double crhs23 = C_3_4*DN_0_1;
const double crhs24 = C_4_5*DN_0_2;
const double crhs25 = tau1*(DN_0_1*crhs16 + DN_0_2*crhs10 + DN_1_1*crhs17 + DN_1_2*crhs11 + DN_2_1*crhs18 + DN_2_2*crhs12 + DN_3_1*crhs19 + DN_3_2*crhs13);
const double crhs26 = C_0_5*DN_0_0;
const double crhs27 = tau1*(DN_0_0*crhs16 + DN_0_2*crhs4 + DN_1_0*crhs17 + DN_1_2*crhs5 + DN_2_0*crhs18 + DN_2_2*crhs6 + DN_3_0*crhs19 + DN_3_2*crhs7);
const double crhs28 = acch_1*rho - fgauss_1;
const double crhs29 = crhs3*rho;
const double crhs30 = C_2_4*DN_0_2;
const double crhs31 = C_1_4*DN_0_1;
const double crhs32 = C_3_5*DN_0_0;
const double crhs33 = acch_2*rho - fgauss_2;
const double crhs34 = tau1*(crhs0 + grad_p_0);
const double crhs35 = rho*tau1;
const double crhs36 = crhs35*(crhs28 + grad_p_1);
const double crhs37 = crhs35*(crhs33 + grad_p_2);
const double crhs38 = C_1_3*DN_1_1;
const double crhs39 = C_2_5*DN_1_2;
const double crhs40 = C_0_3*DN_1_0;
const double crhs41 = C_3_4*DN_1_1;
const double crhs42 = C_4_5*DN_1_2;
const double crhs43 = C_0_5*DN_1_0;
const double crhs44 = C_2_4*DN_1_2;
const double crhs45 = C_1_4*DN_1_1;
const double crhs46 = C_3_5*DN_1_0;
const double crhs47 = crhs34*rho;
const double crhs48 = C_1_3*DN_2_1;
const double crhs49 = C_2_5*DN_2_2;
const double crhs50 = C_0_3*DN_2_0;
const double crhs51 = C_3_4*DN_2_1;
const double crhs52 = C_4_5*DN_2_2;
const double crhs53 = C_0_5*DN_2_0;
const double crhs54 = C_2_4*DN_2_2;
const double crhs55 = C_1_4*DN_2_1;
const double crhs56 = C_3_5*DN_2_0;
const double crhs57 = C_1_3*DN_3_1;
const double crhs58 = C_2_5*DN_3_2;
const double crhs59 = C_0_3*DN_3_0;
const double crhs60 = C_3_4*DN_3_1;
const double crhs61 = C_4_5*DN_3_2;
const double crhs62 = C_0_5*DN_3_0;
const double crhs63 = C_2_4*DN_3_2;
const double crhs64 = C_1_4*DN_3_1;
const double crhs65 = C_3_5*DN_3_0;
rhs[0]=DN_0_0*pgauss - DN_0_0*stress_0 - DN_0_1*stress_3 - DN_0_2*stress_5 - N_0*crhs0 + crhs1*crhs3 - crhs14*(C_0_1*DN_0_0 + C_1_5*DN_0_2 + crhs9) - crhs20*(C_0_2*DN_0_0 + C_2_3*DN_0_1 + crhs15) - crhs22*(C_3_3*DN_0_1 + C_3_5*DN_0_2 + crhs21) - crhs25*(C_0_4*DN_0_0 + crhs23 + crhs24) - crhs27*(C_3_5*DN_0_1 + C_5_5*DN_0_2 + crhs26) - crhs8*(C_0_0*DN_0_0 + C_0_3*DN_0_1 + C_0_5*DN_0_2);
rhs[1]=-DN_0_0*stress_3 + DN_0_1*crhs29 + DN_0_1*pgauss - DN_0_1*stress_1 - DN_0_2*stress_4 - N_0*crhs28 - crhs14*(C_1_1*DN_0_1 + C_1_3*DN_0_0 + C_1_4*DN_0_2) - crhs20*(C_1_2*DN_0_1 + C_2_3*DN_0_0 + crhs30) - crhs22*(C_3_3*DN_0_0 + C_3_4*DN_0_2 + crhs9) - crhs25*(C_3_4*DN_0_0 + C_4_4*DN_0_2 + crhs31) - crhs27*(C_1_5*DN_0_1 + crhs24 + crhs32) - crhs8*(C_0_1*DN_0_1 + C_0_4*DN_0_2 + crhs21);
rhs[2]=-DN_0_0*stress_5 - DN_0_1*stress_4 + DN_0_2*crhs29 + DN_0_2*pgauss - DN_0_2*stress_2 - N_0*crhs33 - crhs14*(C_1_2*DN_0_2 + C_1_5*DN_0_0 + crhs31) - crhs20*(C_2_2*DN_0_2 + C_2_4*DN_0_1 + C_2_5*DN_0_0) - crhs22*(C_2_3*DN_0_2 + crhs23 + crhs32) - crhs25*(C_4_4*DN_0_1 + C_4_5*DN_0_0 + crhs30) - crhs27*(C_4_5*DN_0_1 + C_5_5*DN_0_0 + crhs15) - crhs8*(C_0_2*DN_0_2 + C_0_4*DN_0_1 + crhs26);
rhs[3]=-DN_0_1*crhs36 - DN_0_2*crhs37 - N_0*crhs2 - crhs1*crhs34;
rhs[4]=DN_1_0*crhs29 + DN_1_0*pgauss - DN_1_0*stress_0 - DN_1_1*stress_3 - DN_1_2*stress_5 - N_1*crhs0 - crhs14*(C_0_1*DN_1_0 + C_1_5*DN_1_2 + crhs38) - crhs20*(C_0_2*DN_1_0 + C_2_3*DN_1_1 + crhs39) - crhs22*(C_3_3*DN_1_1 + C_3_5*DN_1_2 + crhs40) - crhs25*(C_0_4*DN_1_0 + crhs41 + crhs42) - crhs27*(C_3_5*DN_1_1 + C_5_5*DN_1_2 + crhs43) - crhs8*(C_0_0*DN_1_0 + C_0_3*DN_1_1 + C_0_5*DN_1_2);
rhs[5]=-DN_1_0*stress_3 + DN_1_1*crhs29 + DN_1_1*pgauss - DN_1_1*stress_1 - DN_1_2*stress_4 - N_1*crhs28 - crhs14*(C_1_1*DN_1_1 + C_1_3*DN_1_0 + C_1_4*DN_1_2) - crhs20*(C_1_2*DN_1_1 + C_2_3*DN_1_0 + crhs44) - crhs22*(C_3_3*DN_1_0 + C_3_4*DN_1_2 + crhs38) - crhs25*(C_3_4*DN_1_0 + C_4_4*DN_1_2 + crhs45) - crhs27*(C_1_5*DN_1_1 + crhs42 + crhs46) - crhs8*(C_0_1*DN_1_1 + C_0_4*DN_1_2 + crhs40);
rhs[6]=-DN_1_0*stress_5 - DN_1_1*stress_4 + DN_1_2*crhs29 + DN_1_2*pgauss - DN_1_2*stress_2 - N_1*crhs33 - crhs14*(C_1_2*DN_1_2 + C_1_5*DN_1_0 + crhs45) - crhs20*(C_2_2*DN_1_2 + C_2_4*DN_1_1 + C_2_5*DN_1_0) - crhs22*(C_2_3*DN_1_2 + crhs41 + crhs46) - crhs25*(C_4_4*DN_1_1 + C_4_5*DN_1_0 + crhs44) - crhs27*(C_4_5*DN_1_1 + C_5_5*DN_1_0 + crhs39) - crhs8*(C_0_2*DN_1_2 + C_0_4*DN_1_1 + crhs43);
rhs[7]=-DN_1_0*crhs47 - DN_1_1*crhs36 - DN_1_2*crhs37 - N_1*crhs2;
rhs[8]=DN_2_0*crhs29 + DN_2_0*pgauss - DN_2_0*stress_0 - DN_2_1*stress_3 - DN_2_2*stress_5 - N_2*crhs0 - crhs14*(C_0_1*DN_2_0 + C_1_5*DN_2_2 + crhs48) - crhs20*(C_0_2*DN_2_0 + C_2_3*DN_2_1 + crhs49) - crhs22*(C_3_3*DN_2_1 + C_3_5*DN_2_2 + crhs50) - crhs25*(C_0_4*DN_2_0 + crhs51 + crhs52) - crhs27*(C_3_5*DN_2_1 + C_5_5*DN_2_2 + crhs53) - crhs8*(C_0_0*DN_2_0 + C_0_3*DN_2_1 + C_0_5*DN_2_2);
rhs[9]=-DN_2_0*stress_3 + DN_2_1*crhs29 + DN_2_1*pgauss - DN_2_1*stress_1 - DN_2_2*stress_4 - N_2*crhs28 - crhs14*(C_1_1*DN_2_1 + C_1_3*DN_2_0 + C_1_4*DN_2_2) - crhs20*(C_1_2*DN_2_1 + C_2_3*DN_2_0 + crhs54) - crhs22*(C_3_3*DN_2_0 + C_3_4*DN_2_2 + crhs48) - crhs25*(C_3_4*DN_2_0 + C_4_4*DN_2_2 + crhs55) - crhs27*(C_1_5*DN_2_1 + crhs52 + crhs56) - crhs8*(C_0_1*DN_2_1 + C_0_4*DN_2_2 + crhs50);
rhs[10]=-DN_2_0*stress_5 - DN_2_1*stress_4 + DN_2_2*crhs29 + DN_2_2*pgauss - DN_2_2*stress_2 - N_2*crhs33 - crhs14*(C_1_2*DN_2_2 + C_1_5*DN_2_0 + crhs55) - crhs20*(C_2_2*DN_2_2 + C_2_4*DN_2_1 + C_2_5*DN_2_0) - crhs22*(C_2_3*DN_2_2 + crhs51 + crhs56) - crhs25*(C_4_4*DN_2_1 + C_4_5*DN_2_0 + crhs54) - crhs27*(C_4_5*DN_2_1 + C_5_5*DN_2_0 + crhs49) - crhs8*(C_0_2*DN_2_2 + C_0_4*DN_2_1 + crhs53);
rhs[11]=-DN_2_0*crhs47 - DN_2_1*crhs36 - DN_2_2*crhs37 - N_2*crhs2;
rhs[12]=DN_3_0*crhs29 + DN_3_0*pgauss - DN_3_0*stress_0 - DN_3_1*stress_3 - DN_3_2*stress_5 - N_3*crhs0 - crhs14*(C_0_1*DN_3_0 + C_1_5*DN_3_2 + crhs57) - crhs20*(C_0_2*DN_3_0 + C_2_3*DN_3_1 + crhs58) - crhs22*(C_3_3*DN_3_1 + C_3_5*DN_3_2 + crhs59) - crhs25*(C_0_4*DN_3_0 + crhs60 + crhs61) - crhs27*(C_3_5*DN_3_1 + C_5_5*DN_3_2 + crhs62) - crhs8*(C_0_0*DN_3_0 + C_0_3*DN_3_1 + C_0_5*DN_3_2);
rhs[13]=-DN_3_0*stress_3 + DN_3_1*crhs29 + DN_3_1*pgauss - DN_3_1*stress_1 - DN_3_2*stress_4 - N_3*crhs28 - crhs14*(C_1_1*DN_3_1 + C_1_3*DN_3_0 + C_1_4*DN_3_2) - crhs20*(C_1_2*DN_3_1 + C_2_3*DN_3_0 + crhs63) - crhs22*(C_3_3*DN_3_0 + C_3_4*DN_3_2 + crhs57) - crhs25*(C_3_4*DN_3_0 + C_4_4*DN_3_2 + crhs64) - crhs27*(C_1_5*DN_3_1 + crhs61 + crhs65) - crhs8*(C_0_1*DN_3_1 + C_0_4*DN_3_2 + crhs59);
rhs[14]=-DN_3_0*stress_5 - DN_3_1*stress_4 + DN_3_2*crhs29 + DN_3_2*pgauss - DN_3_2*stress_2 - N_3*crhs33 - crhs14*(C_1_2*DN_3_2 + C_1_5*DN_3_0 + crhs64) - crhs20*(C_2_2*DN_3_2 + C_2_4*DN_3_1 + C_2_5*DN_3_0) - crhs22*(C_2_3*DN_3_2 + crhs60 + crhs65) - crhs25*(C_4_4*DN_3_1 + C_4_5*DN_3_0 + crhs63) - crhs27*(C_4_5*DN_3_1 + C_5_5*DN_3_0 + crhs58) - crhs8*(C_0_2*DN_3_2 + C_0_4*DN_3_1 + crhs62);
rhs[15]=-DN_3_0*crhs47 - DN_3_1*crhs36 - DN_3_2*crhs37 - N_3*crhs2;

    }


void Stokes3DTwoFluid::ComputeGaussPointEnrichmentContributions(
    BoundedMatrix<double,4,16>& H,
    BoundedMatrix<double,16,4>& V,
    BoundedMatrix<double,4,4>&  Kee,
    array_1d<double,4>& rhs_ee,
    const element_data<4,3>& data,
    const array_1d<double,4>& distances,
    const array_1d<double,4>& Nenr,
    const BoundedMatrix<double,4,4>& DNenr
    )
    {
        const int nnodes = 4;
        const int dim = 3;

        const double rho = inner_prod(data.N, data.rho);
        const double& bdf0 = data.bdf0;
        const double& bdf1 = data.bdf1;
        const double& bdf2 = data.bdf2;

        const BoundedMatrix<double,nnodes,dim>& v = data.v;
        const BoundedMatrix<double,nnodes,dim>& vn = data.vn;
        const BoundedMatrix<double,nnodes,dim>& vnn = data.vnn;
        const BoundedMatrix<double,nnodes,dim>& f = data.f;
        const array_1d<double,nnodes>& p = data.p;

        //get constitutive matrix
        const Matrix& C = data.C;
//         const Vector& stress = data.stress;

        //get shape function values
        const BoundedMatrix<double,nnodes,dim>& DN = data.DN_DX;
        const array_1d<double,nnodes>& N = data.N;

        //compute an equivalent tau by Bitrans*c*Bi
        const double tau_denom = (C(3,3) + C(4,4) + C(5,5))*(pow(DN(0,0), 2) + pow(DN(0,1), 2) + pow(DN(0,2), 2) + pow(DN(1,0), 2) + pow(DN(1,1), 2) + pow(DN(1,2), 2) + pow(DN(2,0), 2) + pow(DN(2,1), 2) + pow(DN(2,2), 2) + pow(DN(3,0), 2) + pow(DN(3,1), 2) + pow(DN(3,2), 2));

        const double tau1 = 1.0/(tau_denom*rho);

        //auxiliary variables used in the calculation of the RHS
        const array_1d<double,dim> fgauss = prod(trans(f), N);
        const array_1d<double,dim> vgauss = prod(trans(v), N);
        const array_1d<double,dim> grad_p = prod(trans(DN), p);
//         const double pgauss = inner_prod(N,p);

        array_1d<double,dim> acch = bdf0*vgauss;
        noalias(acch) += bdf1*prod(trans(vn), N);
        noalias(acch) += bdf2*prod(trans(vnn), N);

        array_1d<double,4> penr = ZeroVector(4); //penriched is considered to be zero as we do not want to store it

        const double DN_0_0 = DN(0,0); const double DN_0_1 = DN(0,1); const double DN_0_2 = DN(0,2);
        const double DN_1_0 = DN(1,0); const double DN_1_1 = DN(1,1); const double DN_1_2 = DN(1,2);
        const double DN_2_0 = DN(2,0); const double DN_2_1 = DN(2,1); const double DN_2_2 = DN(2,2);
        const double DN_3_0 = DN(3,0); const double DN_3_1 = DN(3,1); const double DN_3_2 = DN(3,2);

        const double Nenr_0 = Nenr[0];
        const double Nenr_1 = Nenr[1];
        const double Nenr_2 = Nenr[2];
        const double Nenr_3 = Nenr[3];

        const double N_0 = N[0];
        const double N_1 = N[1];
        const double N_2 = N[2];
        const double N_3 = N[3];

        const double p_0 = p[0];
        const double p_1 = p[1];
        const double p_2 = p[2];
        const double p_3 = p[3];
        
        const double penr_0 = penr[0];
        const double penr_1 = penr[1];
        const double penr_2 = penr[2];
        const double penr_3 = penr[3];

        const double DNenr_0_0 = DNenr(0,0); const double DNenr_0_1 = DNenr(0,1); const double DNenr_0_2 = DNenr(0,2);
        const double DNenr_1_0 = DNenr(1,0); const double DNenr_1_1 = DNenr(1,1); const double DNenr_1_2 = DNenr(1,2);
        const double DNenr_2_0 = DNenr(2,0); const double DNenr_2_1 = DNenr(2,1); const double DNenr_2_2 = DNenr(2,2);
        const double DNenr_3_0 = DNenr(3,0); const double DNenr_3_1 = DNenr(3,1); const double DNenr_3_2 = DNenr(3,2);

        const double f_0_0 = f(0,0); const double f_0_1 = f(0,1); const double f_0_2 = f(0,2);
        const double f_1_0 = f(1,0); const double f_1_1 = f(1,1); const double f_1_2 = f(1,2);
        const double f_2_0 = f(2,0); const double f_2_1 = f(2,1); const double f_2_2 = f(2,2);
        const double f_3_0 = f(3,0); const double f_3_1 = f(3,1); const double f_3_2 = f(3,2);

        const double v_0_0 = v(0,0); const double v_0_1 = v(0,1); const double v_0_2 = v(0,2); 
        const double v_1_0 = v(1,0); const double v_1_1 = v(1,1); const double v_1_2 = v(1,2);
        const double v_2_0 = v(2,0); const double v_2_1 = v(2,1); const double v_2_2 = v(2,2);
        const double v_3_0 = v(3,0); const double v_3_1 = v(3,1); const double v_3_2 = v(3,2);

        const double vn_0_0 = vn(0,0); const double vn_0_1 = vn(0,1); const double vn_0_2 = vn(0,2); 
        const double vn_1_0 = vn(1,0); const double vn_1_1 = vn(1,1); const double vn_1_2 = vn(1,2);
        const double vn_2_0 = vn(2,0); const double vn_2_1 = vn(2,1); const double vn_2_2 = vn(2,2);
        const double vn_3_0 = vn(3,0); const double vn_3_1 = vn(3,1); const double vn_3_2 = vn(3,2);

        const double vnn_0_0 = vnn(0,0); const double vnn_0_1 = vnn(0,1); const double vnn_0_2 = vnn(0,2); 
        const double vnn_1_0 = vnn(1,0); const double vnn_1_1 = vnn(1,1); const double vnn_1_2 = vnn(1,2);
        const double vnn_2_0 = vnn(2,0); const double vnn_2_1 = vnn(2,1); const double vnn_2_2 = vnn(2,2);
        const double vnn_3_0 = vnn(3,0); const double vnn_3_1 = vnn(3,1); const double vnn_3_2 = vnn(3,2);

        const double cV0 = rho*tau1;
V(0,0)=-DN_0_0*Nenr_0;
V(0,1)=-DN_0_0*Nenr_1;
V(0,2)=-DN_0_0*Nenr_2;
V(0,3)=-DN_0_0*Nenr_3;
V(1,0)=-DN_0_1*Nenr_0;
V(1,1)=-DN_0_1*Nenr_1;
V(1,2)=-DN_0_1*Nenr_2;
V(1,3)=-DN_0_1*Nenr_3;
V(2,0)=-DN_0_2*Nenr_0;
V(2,1)=-DN_0_2*Nenr_1;
V(2,2)=-DN_0_2*Nenr_2;
V(2,3)=-DN_0_2*Nenr_3;
V(3,0)=cV0*(DN_0_0*DNenr_0_0 + DN_0_1*DNenr_0_1 + DN_0_2*DNenr_0_2);
V(3,1)=cV0*(DN_0_0*DNenr_1_0 + DN_0_1*DNenr_1_1 + DN_0_2*DNenr_1_2);
V(3,2)=cV0*(DN_0_0*DNenr_2_0 + DN_0_1*DNenr_2_1 + DN_0_2*DNenr_2_2);
V(3,3)=cV0*(DN_0_0*DNenr_3_0 + DN_0_1*DNenr_3_1 + DN_0_2*DNenr_3_2);
V(4,0)=-DN_1_0*Nenr_0;
V(4,1)=-DN_1_0*Nenr_1;
V(4,2)=-DN_1_0*Nenr_2;
V(4,3)=-DN_1_0*Nenr_3;
V(5,0)=-DN_1_1*Nenr_0;
V(5,1)=-DN_1_1*Nenr_1;
V(5,2)=-DN_1_1*Nenr_2;
V(5,3)=-DN_1_1*Nenr_3;
V(6,0)=-DN_1_2*Nenr_0;
V(6,1)=-DN_1_2*Nenr_1;
V(6,2)=-DN_1_2*Nenr_2;
V(6,3)=-DN_1_2*Nenr_3;
V(7,0)=cV0*(DN_1_0*DNenr_0_0 + DN_1_1*DNenr_0_1 + DN_1_2*DNenr_0_2);
V(7,1)=cV0*(DN_1_0*DNenr_1_0 + DN_1_1*DNenr_1_1 + DN_1_2*DNenr_1_2);
V(7,2)=cV0*(DN_1_0*DNenr_2_0 + DN_1_1*DNenr_2_1 + DN_1_2*DNenr_2_2);
V(7,3)=cV0*(DN_1_0*DNenr_3_0 + DN_1_1*DNenr_3_1 + DN_1_2*DNenr_3_2);
V(8,0)=-DN_2_0*Nenr_0;
V(8,1)=-DN_2_0*Nenr_1;
V(8,2)=-DN_2_0*Nenr_2;
V(8,3)=-DN_2_0*Nenr_3;
V(9,0)=-DN_2_1*Nenr_0;
V(9,1)=-DN_2_1*Nenr_1;
V(9,2)=-DN_2_1*Nenr_2;
V(9,3)=-DN_2_1*Nenr_3;
V(10,0)=-DN_2_2*Nenr_0;
V(10,1)=-DN_2_2*Nenr_1;
V(10,2)=-DN_2_2*Nenr_2;
V(10,3)=-DN_2_2*Nenr_3;
V(11,0)=cV0*(DN_2_0*DNenr_0_0 + DN_2_1*DNenr_0_1 + DN_2_2*DNenr_0_2);
V(11,1)=cV0*(DN_2_0*DNenr_1_0 + DN_2_1*DNenr_1_1 + DN_2_2*DNenr_1_2);
V(11,2)=cV0*(DN_2_0*DNenr_2_0 + DN_2_1*DNenr_2_1 + DN_2_2*DNenr_2_2);
V(11,3)=cV0*(DN_2_0*DNenr_3_0 + DN_2_1*DNenr_3_1 + DN_2_2*DNenr_3_2);
V(12,0)=-DN_3_0*Nenr_0;
V(12,1)=-DN_3_0*Nenr_1;
V(12,2)=-DN_3_0*Nenr_2;
V(12,3)=-DN_3_0*Nenr_3;
V(13,0)=-DN_3_1*Nenr_0;
V(13,1)=-DN_3_1*Nenr_1;
V(13,2)=-DN_3_1*Nenr_2;
V(13,3)=-DN_3_1*Nenr_3;
V(14,0)=-DN_3_2*Nenr_0;
V(14,1)=-DN_3_2*Nenr_1;
V(14,2)=-DN_3_2*Nenr_2;
V(14,3)=-DN_3_2*Nenr_3;
V(15,0)=cV0*(DN_3_0*DNenr_0_0 + DN_3_1*DNenr_0_1 + DN_3_2*DNenr_0_2);
V(15,1)=cV0*(DN_3_0*DNenr_1_0 + DN_3_1*DNenr_1_1 + DN_3_2*DNenr_1_2);
V(15,2)=cV0*(DN_3_0*DNenr_2_0 + DN_3_1*DNenr_2_1 + DN_3_2*DNenr_2_2);
V(15,3)=cV0*(DN_3_0*DNenr_3_0 + DN_3_1*DNenr_3_1 + DN_3_2*DNenr_3_2);


        const double cH0 = bdf0*pow(rho, 2)*tau1;
const double cH1 = N_0*cH0;
const double cH2 = rho*tau1;
const double cH3 = N_1*cH0;
const double cH4 = N_2*cH0;
const double cH5 = N_3*cH0;
H(0,0)=DNenr_0_0*cH1;
H(0,1)=DNenr_0_1*cH1;
H(0,2)=DNenr_0_2*cH1;
H(0,3)=cH2*(DN_0_0*DNenr_0_0 + DN_0_1*DNenr_0_1 + DN_0_2*DNenr_0_2);
H(0,4)=DNenr_0_0*cH3;
H(0,5)=DNenr_0_1*cH3;
H(0,6)=DNenr_0_2*cH3;
H(0,7)=cH2*(DN_1_0*DNenr_0_0 + DN_1_1*DNenr_0_1 + DN_1_2*DNenr_0_2);
H(0,8)=DNenr_0_0*cH4;
H(0,9)=DNenr_0_1*cH4;
H(0,10)=DNenr_0_2*cH4;
H(0,11)=cH2*(DN_2_0*DNenr_0_0 + DN_2_1*DNenr_0_1 + DN_2_2*DNenr_0_2);
H(0,12)=DNenr_0_0*cH5;
H(0,13)=DNenr_0_1*cH5;
H(0,14)=DNenr_0_2*cH5;
H(0,15)=cH2*(DN_3_0*DNenr_0_0 + DN_3_1*DNenr_0_1 + DN_3_2*DNenr_0_2);
H(1,0)=DNenr_1_0*cH1;
H(1,1)=DNenr_1_1*cH1;
H(1,2)=DNenr_1_2*cH1;
H(1,3)=cH2*(DN_0_0*DNenr_1_0 + DN_0_1*DNenr_1_1 + DN_0_2*DNenr_1_2);
H(1,4)=DNenr_1_0*cH3;
H(1,5)=DNenr_1_1*cH3;
H(1,6)=DNenr_1_2*cH3;
H(1,7)=cH2*(DN_1_0*DNenr_1_0 + DN_1_1*DNenr_1_1 + DN_1_2*DNenr_1_2);
H(1,8)=DNenr_1_0*cH4;
H(1,9)=DNenr_1_1*cH4;
H(1,10)=DNenr_1_2*cH4;
H(1,11)=cH2*(DN_2_0*DNenr_1_0 + DN_2_1*DNenr_1_1 + DN_2_2*DNenr_1_2);
H(1,12)=DNenr_1_0*cH5;
H(1,13)=DNenr_1_1*cH5;
H(1,14)=DNenr_1_2*cH5;
H(1,15)=cH2*(DN_3_0*DNenr_1_0 + DN_3_1*DNenr_1_1 + DN_3_2*DNenr_1_2);
H(2,0)=DNenr_2_0*cH1;
H(2,1)=DNenr_2_1*cH1;
H(2,2)=DNenr_2_2*cH1;
H(2,3)=cH2*(DN_0_0*DNenr_2_0 + DN_0_1*DNenr_2_1 + DN_0_2*DNenr_2_2);
H(2,4)=DNenr_2_0*cH3;
H(2,5)=DNenr_2_1*cH3;
H(2,6)=DNenr_2_2*cH3;
H(2,7)=cH2*(DN_1_0*DNenr_2_0 + DN_1_1*DNenr_2_1 + DN_1_2*DNenr_2_2);
H(2,8)=DNenr_2_0*cH4;
H(2,9)=DNenr_2_1*cH4;
H(2,10)=DNenr_2_2*cH4;
H(2,11)=cH2*(DN_2_0*DNenr_2_0 + DN_2_1*DNenr_2_1 + DN_2_2*DNenr_2_2);
H(2,12)=DNenr_2_0*cH5;
H(2,13)=DNenr_2_1*cH5;
H(2,14)=DNenr_2_2*cH5;
H(2,15)=cH2*(DN_3_0*DNenr_2_0 + DN_3_1*DNenr_2_1 + DN_3_2*DNenr_2_2);
H(3,0)=DNenr_3_0*cH1;
H(3,1)=DNenr_3_1*cH1;
H(3,2)=DNenr_3_2*cH1;
H(3,3)=cH2*(DN_0_0*DNenr_3_0 + DN_0_1*DNenr_3_1 + DN_0_2*DNenr_3_2);
H(3,4)=DNenr_3_0*cH3;
H(3,5)=DNenr_3_1*cH3;
H(3,6)=DNenr_3_2*cH3;
H(3,7)=cH2*(DN_1_0*DNenr_3_0 + DN_1_1*DNenr_3_1 + DN_1_2*DNenr_3_2);
H(3,8)=DNenr_3_0*cH4;
H(3,9)=DNenr_3_1*cH4;
H(3,10)=DNenr_3_2*cH4;
H(3,11)=cH2*(DN_2_0*DNenr_3_0 + DN_2_1*DNenr_3_1 + DN_2_2*DNenr_3_2);
H(3,12)=DNenr_3_0*cH5;
H(3,13)=DNenr_3_1*cH5;
H(3,14)=DNenr_3_2*cH5;
H(3,15)=cH2*(DN_3_0*DNenr_3_0 + DN_3_1*DNenr_3_1 + DN_3_2*DNenr_3_2);


        const double cKee0 = rho*tau1;
const double cKee1 = cKee0*(DNenr_0_0*DNenr_1_0 + DNenr_0_1*DNenr_1_1 + DNenr_0_2*DNenr_1_2);
const double cKee2 = cKee0*(DNenr_0_0*DNenr_2_0 + DNenr_0_1*DNenr_2_1 + DNenr_0_2*DNenr_2_2);
const double cKee3 = cKee0*(DNenr_0_0*DNenr_3_0 + DNenr_0_1*DNenr_3_1 + DNenr_0_2*DNenr_3_2);
const double cKee4 = cKee0*(DNenr_1_0*DNenr_2_0 + DNenr_1_1*DNenr_2_1 + DNenr_1_2*DNenr_2_2);
const double cKee5 = cKee0*(DNenr_1_0*DNenr_3_0 + DNenr_1_1*DNenr_3_1 + DNenr_1_2*DNenr_3_2);
const double cKee6 = cKee0*(DNenr_2_0*DNenr_3_0 + DNenr_2_1*DNenr_3_1 + DNenr_2_2*DNenr_3_2);
Kee(0,0)=cKee0*(pow(DNenr_0_0, 2) + pow(DNenr_0_1, 2) + pow(DNenr_0_2, 2));
Kee(0,1)=cKee1;
Kee(0,2)=cKee2;
Kee(0,3)=cKee3;
Kee(1,0)=cKee1;
Kee(1,1)=cKee0*(pow(DNenr_1_0, 2) + pow(DNenr_1_1, 2) + pow(DNenr_1_2, 2));
Kee(1,2)=cKee4;
Kee(1,3)=cKee5;
Kee(2,0)=cKee2;
Kee(2,1)=cKee4;
Kee(2,2)=cKee0*(pow(DNenr_2_0, 2) + pow(DNenr_2_1, 2) + pow(DNenr_2_2, 2));
Kee(2,3)=cKee6;
Kee(3,0)=cKee3;
Kee(3,1)=cKee5;
Kee(3,2)=cKee6;
Kee(3,3)=cKee0*(pow(DNenr_3_0, 2) + pow(DNenr_3_1, 2) + pow(DNenr_3_2, 2));


        const double crhs_ee0 = DN_0_0*p_0 + DN_1_0*p_1 + DN_2_0*p_2 + DN_3_0*p_3 + DNenr_0_0*penr_0 + DNenr_1_0*penr_1 + DNenr_2_0*penr_2 + DNenr_3_0*penr_3 - N_0*f_0_0 - N_1*f_1_0 - N_2*f_2_0 - N_3*f_3_0 + rho*(N_0*(bdf0*v_0_0 + bdf1*vn_0_0 + bdf2*vnn_0_0) + N_1*(bdf0*v_1_0 + bdf1*vn_1_0 + bdf2*vnn_1_0) + N_2*(bdf0*v_2_0 + bdf1*vn_2_0 + bdf2*vnn_2_0) + N_3*(bdf0*v_3_0 + bdf1*vn_3_0 + bdf2*vnn_3_0));
const double crhs_ee1 = DN_0_1*p_0 + DN_1_1*p_1 + DN_2_1*p_2 + DN_3_1*p_3 + DNenr_0_1*penr_0 + DNenr_1_1*penr_1 + DNenr_2_1*penr_2 + DNenr_3_1*penr_3 - N_0*f_0_1 - N_1*f_1_1 - N_2*f_2_1 - N_3*f_3_1 + rho*(N_0*(bdf0*v_0_1 + bdf1*vn_0_1 + bdf2*vnn_0_1) + N_1*(bdf0*v_1_1 + bdf1*vn_1_1 + bdf2*vnn_1_1) + N_2*(bdf0*v_2_1 + bdf1*vn_2_1 + bdf2*vnn_2_1) + N_3*(bdf0*v_3_1 + bdf1*vn_3_1 + bdf2*vnn_3_1));
const double crhs_ee2 = DN_0_2*p_0 + DN_1_2*p_1 + DN_2_2*p_2 + DN_3_2*p_3 + DNenr_0_2*penr_0 + DNenr_1_2*penr_1 + DNenr_2_2*penr_2 + DNenr_3_2*penr_3 - N_0*f_0_2 - N_1*f_1_2 - N_2*f_2_2 - N_3*f_3_2 + rho*(N_0*(bdf0*v_0_2 + bdf1*vn_0_2 + bdf2*vnn_0_2) + N_1*(bdf0*v_1_2 + bdf1*vn_1_2 + bdf2*vnn_1_2) + N_2*(bdf0*v_2_2 + bdf1*vn_2_2 + bdf2*vnn_2_2) + N_3*(bdf0*v_3_2 + bdf1*vn_3_2 + bdf2*vnn_3_2));
const double crhs_ee3 = rho*tau1;
rhs_ee[0]=-crhs_ee3*(DNenr_0_0*crhs_ee0 + DNenr_0_1*crhs_ee1 + DNenr_0_2*crhs_ee2);
rhs_ee[1]=-crhs_ee3*(DNenr_1_0*crhs_ee0 + DNenr_1_1*crhs_ee1 + DNenr_1_2*crhs_ee2);
rhs_ee[2]=-crhs_ee3*(DNenr_2_0*crhs_ee0 + DNenr_2_1*crhs_ee1 + DNenr_2_2*crhs_ee2);
rhs_ee[3]=-crhs_ee3*(DNenr_3_0*crhs_ee0 + DNenr_3_1*crhs_ee1 + DNenr_3_2*crhs_ee2);



    }

}